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Plan of the thesis 



During the last five years, I have been learning some of the theoretical tools of 
Many-Body physics, and I have tried to apply them to the fascinating world 
of quantum systems. Thus, the studies presented in this thesis constitute an 
attempt to show the usefulness of these techniques for the understanding of a 
variety of systems ranging from ultracold gases (both fermionic and bosonic) 
to the more familiar helium. The range of techniques that we have used 
is comparable to the number of systems under study: from the simple, but 
always illustrative mean-field approach, to the powerful Green's functions' 
formalism, and from Monte Carlo methods to the Density Functional theory. 

I have made an effort to present the various works in an understand- 
able way. I would like to thank Prof. Artur Polls, my advisor, and Dr. 
Lloreng Brualla, Dr. Muntsa Guilleumas, Prof. Jesus Navarro and Dr. Ar- 
men Sedrakian for carefully reading various parts of the manuscript, raising 
interesting questions and suggesting improvements. In any case, of course, 
any error or misunderstanding should be attributed only to me. I'll be glad 
to receive questions or comments on any point related to this work.* 

The plan of the thesis is as follows. In a first part, we study ultracold 
gases of alkali atoms. The first three chapters are devoted to the evaluation 
of the possibilities of having and detecting a superfiuid system of fermionic 
atoms. The first chapter introduces the idea of pairing in a Green's func- 
tions' formulation of the theory of superconductivity by Bardeen, Cooper and 
Schrieffer. Chapter [21 presents the application of this theory to the case of a 
two-component system where the two species may have different densities, a 
fact that has important consequences for the prospects of superfiuidity. This 
work was performed in collaboration with Dr. Hans- Josef Schulze who, at 
that time, was a post-doc in Barcelona. Many of the results expounded here 
have appeared in the article |MPSfll] . 

Chapter 01 generalizes the structure of the BCS ground state to the cases 
of finite-momentum Cooper pairs (the so-called loff or fflo phase) and 
deformed Fermi surfaces (dfs). This last state was first proposed by Prof. 



*You can contact me at jordimp2003@yahoo.es 
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Herbert Miither and Dr. Armen Sedrakian from Universitat Tubingen, with 
whom we have worked together and who I would like to thank for their warm 
hospitality during my visits to Tiibingen. As a result of this work, we have 
recently published the paper |SMPMfl5] . 

Chapter |3| introduces for the first time bosons into this thesis. They 
are used to cool a one-component system of fermions, and their capacity 
to help the fermions to pair is studied, with emphasis in the case of a two- 
dimensional configuration. This research, which appeared in |MPBSn4b| and 
m |MPRS04a| . was performed again with Dr. Schulze, then in Catania, and 
with Prof. Marcello Baldo. I want to thank both of them, and the rest of 
the INFN group for their friendship and hospitality during my visit in 2003. 

In chapter Elwe enter boldly into the world of bosonic systems and analyze 
the dynamical evolution of a set of bosons whose spin degree of freedom 
is free to evolve inside an optical trap, only restricted by the conservation 
of magnetization. An attempt to incorporate temperature effects is also 
performed. Thanks to this project I could enjoy for a couple of months the 
stimulating atmosphere in the group of Prof. Maciej Lewenstein and Prof. 
Anna Sanpera at Universitat Hannover, just before they moved to Catalonia. 
This work, developed in close collaboration with the experimental group of 
Prof. K. Sengstock and Dr. K. Bongs from Universitat Hamburg, has been 
submitted for publication to Physical Review A |MGS"'"05| . 

At the conceptual half of the thesis (p. IHH), we leave ultracold gases 
and start the study of helium-4 in two dimensions, which constitutes the 
second part of the work. First, in chapter we present a Monte Carlo 
study of the ground-state structure and energetics of "^He puddles, with an 
estimation of the line tension. The results obtained here are incorporated in 
chapter m to build a Density Functional appropiate to study two-dimensional 
helium systems. This functional is then used to analyze a variety of such 
systems, and a comparison with previous studies in two and three dimensions 
is presented. Both these works have been developed in collaboration with 
Prof. Jesus Navarro, from CSIC-Universitat de Valencia, and Dr. Antonio 
Sarsa, who is now at Universidad de Cordoba. Some of the results presented 
here have appeared in |SMPN03| and |MSNP05| . 

The conclusions of the thesis are finally summarized. 



Book I 
Ultracold gases 



Chapter 1 

The pairing solution 



— Bien parece — respondio don Quijote — que no estas cursado en 
esto de las aventuras: ellos son gigantes; y si tienes miedo, quitate 
de ahi, y ponte en oracion en el espacio que yo voy a entrar con 
ellos en fiera y desigual batalla. (...) jNon fuyades, cobardes y viles 
criaturas, que un solo caballero es el que os aconnete! 

Miguel de Cervantes, Don Quijote de la Mancha (I, 8) 

1.1 Historical background: before BCS theory 

The BCS theory of superconductivity |BCS571 ISch88j has been one of the 
most successful contributions to physics in the 20**^ century, because it was 
the first microscopic theory of a truly macroscopic quantum phenomenon. 
It explains the mechanisms behind dissipation-free electric transport in a 
number of materials. Indeed, since its formulation, this theory has been 
applied to such different systems as solid metals and alloys |Sch88[ ldG66j . 
atomic nuclei and neutron stars |BMP58jlDHfl3j . elementary particle physics 
(see, e. g., \Cm4\ ) and superfiuid ^He |()RL72[ [LejTSl . 

The theory developed by Bardeen, Cooper and Schrieffer appeared as a 
final step in the theoretical understanding of a long series of experimental 
discoveries. In 1911 Heike Kamerling Onnes |Kamll| observed that mercury 
below 4.2 K looses its electrical resistance and enters 'a new state, which, 
owing to its particular electrical properties, can be called the state of super- 
conductivity' |Kam67] . In 1933 Meissner and Ochsenfeld |M033] observed 
that a superconductor has also notable magnetic properties, namely it is a 
perfect diamagnet: a (small) applied magnetic field vanishes in the interior 
of a bulk superconductor. Another important contribution was the discovery 
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in 1950 of the isotope effect («. e., the dependence of the critical tempera- 
ture on the mass of the ions of a superconducting solid) by Maxwell [M axSOj 
and Reynolds et al. |R,SWN5fl] . Finally, let us mention the experimental 
determination of the quantization of the magnetic flux traversing a multiply- 
connected superconductor by Deaver and Fairbank |DF61| and Doll and 
Nabauer |DN61| following an initial prediction of the London model. 

There were many theoretical attemps to explain the experimental results 
prior to the theory of Bardeen, Cooper and Schrieffer. In a first phenomeno- 
logical theory, Gorter and Casimir |GC34| formulated a two-fiuid model (sim- 
ilar to that of Tisza |Tis38[ ITis4n| and Landau |Lan41| for "^He) where elec- 
trons can be either in a superconducting state or in the normal state. At 
T = all electrons are in the superfiuid, while the fraction of normal electrons 
grows with temperature and finally equals one at the critical temperatre T^. 
The fraction of normal electrons was calculated by minimazing a free energy 
interpolated between the limits of a superfiuid system at T = and a nor- 
mal one at T = Tc. This model, however, had a number of artificial points 
(such as the way to construct the free energy), just intended to reproduce 
the experimental results known at the moment, but without a microscopic 
physical motivation. 

A more successful theory was that of F. London and H. London |LL35I 
ILon35| based on the electromagnetic phenomenology of supercondutors. It 
explained the Meissner-Ochsenfeld effect [MQ33J. introduced a 'penetration 
length' A of magnetic fields into superconductors and also predicted the quan- 
tization of magnetic fiux in a multiply-connected superfiuid. F. London's 
work also speculated on the possible existence of a gap in the excitation 
spectrum of a superfiuid, a crucial ingredient in BCS theory. 

In 1950 Ginzburg and Landau |LG50] formulated a non-local theory of 
superconductivity by generalizing the ideas of the London brothers: they 
introduced a spatially-varying 'order parameter' closely related to Londons' 
superfiuid density. This theory was valid only near T^, but nevertheless it 
became of great interest with the discovery of type II superconductors. It is 
worth noticing that Ginzburg and Landau's theory can be also derived from 
the BCS microscopic theory |(Tor59| . 

In 1956, Cooper |Coo56| studied the problem of an interacting pair of 
fermions above a frozen Fermi sea and showed that, for an attractive in- 
teraction (no matter how weak it was), the pair would be bound, with an 
exponentially small binding energy in the limit of small coupling. These 
bound pairs are called 'Cooper pairs', and are a central ingredient of the 
BCS theory of superconductivity: Assuming that all the electrons in a metal 
are forming such pairs one can show that the total energy of the system is 
lowered with respect to that of the normal system. 



1.2 The Bardeen-Cooper-SchriefFer theory 
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1.2 The Bardeen-Cooper-Schrieffer theory 

In this section, we introduce the microscopic theory of superconductivity of 
Bardeen, Cooper and Schrieffer (know as 'bcs theory') in a way that will 
later allow us to apply it to a variety of atomic gases as we shall do in the 
following chapters. We will start in Sect. 11.2.11 by introducing the idea of 
Cooper pairs and show how the presence of a Fermi sea allows an interacting 
pair of fermions to bind itself at arbitrary small attraction. This many- 
body effect (in the sense that the pair would not be bound in the absence 
of the Fermi sea unless the attraction was greater than a minimum value) is 
the essence of BCS theory, which we present in a formal way in Sect. 11.2.21 
Finally, we present some basic results of the theory in Sect. 11.2.31 

1.2.1 Premier: Cooper pairs 

Consider a pair of fermions of mass m in homogeneous space above a 'frozen' 
Fermi sea of non-interacting particles with Fermi energy ep = h'^kp/2m. Let 
us assume that these two particles are distinguishable (for example, they can 
be electrons with different spins, atoms of the same chemical element with 
different hyperfine spins, neutrons and protons, or quarks of different colors), 
and that they interact through a two-body spin-independent potential. The 
presence of the Fermi sea forbids the particles from occupying the energy 
levels below ep (see Fig. ll.ljl . Therefore, we can take as their zero-energy 
level. This problem was originally solved by Cooper |Coo56j . so it is usually 
called the 'Cooper problem' and the bound pairs are known as 'Cooper pairs'. 
The solution we describe is based on f SchS&J . 

The wave function of a pair with momentum q in the center-of-mass 
system (c.m.s.) can be written 

ip{xi, X2) = (t)q{r) exp{iqr ■ R} , 

where we have defined the centre-of-mass coordinate R = {xi + X2)/2 and 
the relative coordinate r = xi — X2. Taking for simplicity q = 0, we have 

^ = 0(r) = = E . (1.1) 

k>kp k>kp 

Here we have applied Pauli's exclusion principle to restrict the summations 
over momenta to the region outside the filled Fermi sea. This wave function 
can be interpreted as a superposition of configurations with single-particle 
wave functions with momenta \k) and \ — k). 
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Figure 1.1: The Cooper problem: two fermions in states \k) and | — k) 
interact above a frozen Fermi sea (filled sphere). The fermions form a pair 
(dashed line) with center-of-mass momentum q = and can scatter into 
states \k') and | — k') for |A;'| > kp. 




The corresponding Schrodinger equation can be written 

{E - Ho)il) = Vil) {E- 2eu)ak = ^ Vkk'ak' , 

k'>kF 

Vkk' ■■= {k,-k\V\k',-k') , 



(1.2) 



where we have introduced the single-particle spectrum corresponding to 
Ho = -h^V^/{2m). 

If we take the simple case of a constant, attractive interaction constrained 
to a shell of width Q above the Fermi surface, that is. 



^ ^ f -|A| Bp < Cfc, Cfc/ < £f + n 

\ otherwise 
then the lowest eigenenergy is 

E = 2sf 



29, 



exp 



1 



2eF - 2Vte "("f^w , 



(1.3) 



where y{kF) — mkF/{2Ti'^h^) is the density of single-particle states for a 
single spin projection at the Fermi surface. The last result is valid in the 
weak-coupling approximation, i/(A;i?)|A| -C 1. 
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Therefore, the pair will be bound {E < 2eF) no matter how small |A| is. 
Note that this is a truly many-body effect as, in the absence of the underlying, 
filled Fermi sea, the pair would not bind for too weak an interaction. Fur- 
thermore, the non-analytic behavior in the weak-coupling limit shows that 
this result could not be found by perturbation theory, as it corresponds to a 
deep modification of the ground state as compared to the ideal Fermi gas. 

1.2.2 Formalism for the theory 

Let us study now a more realistic problem in which we have an interacting 
system of fermions with an internal degree of freedom, that we indicate by a 
(pseudo)spin variable s =T, i- This problem can be modeled by the following 
second-quantized Hamiltonian 

+ ^jd^x j d^x' V{x - x')i^\{x)'4)]{x')i)^{x')i)^{x) (1.4) 

Here, ips{x) (^J(a;)) is the operator that destroys (creates) a particle of 
(pseudo)spin s at position x, and V{x — x') is the interaction potential 
between the fermions. We will be interested in low-density, low-temperature 
systems where the most important contribution to interactions is the s-wave 
one. 

We look for a solution to this problem using the tools of quantum field 
theory, in particular, the Green's functions' formalism. We start linearizing 
the interaction and transforming it into a one-body operator by means of a 
Hartree-Fock-like approximation, 

+ (^t(^)V'-.(^'))^L(=«')4(a^)} • (1.5) 

The first term inside the curly brackets is the direct (or Hartree) term, while 
the second one is the exchange (or Fock) one. In the next step, we consider 
a new term in the linearization process which accounts for the possibility, 
outlined in the previous section, that two fermions with opposite spins form 
a bound pair (cf. Sec. 51 in |FW71| or Sec. 7-2 in |Sch88| ). Even though H 
commutes with the number operator N, this pairing possibility is more easily 
implemented if we work in the grand-canonical ensemble. Physically, we can 
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understand this in the following way: for a large number of particles N in the 
system, the ground state energies of the system with N and A^±l particles are 
nearly degenerate if we substract from them the chemical potential. Then, it 
is reasonable that a wave function in Fock space, with only the average value 
N = (N) fixed, is more flexible from a variational point of view than one with 
a fixed number of particles. In this way, the function can be closer to the 
true ground state of the system immersed in the 'bath' of condensed pairs. 
In fact, it is customary to assume that the Hartree-Fock potential above is 
the same in the normal and superfiuid states, so that we can disregard it 
from now on as its inclusion would not affect the comparison between both 
states. Therefore, we will work with the following effective (grand-canonical) 
Hamiltonian 



X [(^j(x)7A|(a^'))^T(^')^iN +^{(a^)^}(a^') {m^')^^: 

The angular brackets (■ ■ ■) denote a thermal average with K: 



x')x 



(1.6) 



(■■■): = 



Tr [e 



-I3K 



] 



Tr [e-P^] 



where ks stands for Boltzmann's constant. 

Introducing now the imaginary-time variable r G [0, h[3] C 
pendixEj) and the field operators in the Heisenberg picture, 

4(a;r) := e^"/^4(a;)e-^"/^ 
one can show that they satisfy these equations of motion: 



(see Ap- 

(1.7a) 
(1.7b) 



.<9^/'t(a;r) 



■v^ - /it i^^{xT) 



dr \ 2m 

+ J d'yV{y-x) (^t(a^)V'i(2/))V'I( 

',t 



[1.8a) 



h- 



dip^ {xT 



8r ' V 2m^ 

d?yV{x-y) U\{x)i)\{y)) i)^{yT) . 



(1.8b) 
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The solution of the problem is more naturally found by introducing the 
temperature Green's function (or normal propagator) G corresponding to 
each species s =t, |, 



GJxt: x't') :- 



(1.9) 



where is the imaginary-time ordering operator |FW71| . The equation of 
motion for is easily seen to be 



h- 



dG^{xT-x'T') 



dr 



h6{T - t')6{x -x')-[ -— - /it G^{xt; x't') 



2m 



+ I d'yV{y-x){ij^{x)^^{y)){-T, ^l{yT)^Ux'T') ) . (1.10) 



Defining the anomalous propagators F and by 



F{xt; x't) 
F\xt]x't') 



-Tr 



i>\{xT)i)\{x'T') 



(1.11a) 
(1.11b) 



we can rewrite Eq. (jl.lfljl as 



d_ 



2m 



G^{xt; x't' 



+ i jd'yA{x,y)F\yT;x'T') = h6{T-T')6{x~x'), (1.12) 



where we defined the two-point gap function 



A(a;, y) := V{y - x) U^{x)t[j^{y) ) V = ~V{y - x)F{xt^, yT)V . (1.13) 



Here we wrote explicitly the volume V of the system to make clear that 
/S.{x,y) has dimensions of energy. Note that all propagators in space coor- 
dinates have dimensions of density. Moreover, note that F and are not 
operators, neither is one the Hermitian adjoint of the other (as the f symbol 
might seem to indicate): they are just c-functions, as can be seen from their 
definitions as expectations values. 
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In a similar way, we can find the equations of motion for Gj, and F: 



h 



d_ 



2m 



Gi{x't'-xt) 



+ ^ j d'y A*{y,x)Fiyr;x'T') = h5ir - r')5{x - x') , (1.14) 



h 



d_ 



2m 
2m 



F\xT] x't') 



j '^A\v,x)G^{yT-x'r' 



V^-Z^t) F{xr-x'r') = j '^A{x,y)G^{x'r'-yr) 



(1.15) 



[1.16) 



Then, in the absence of external fields, translational symmetry implies 
that a Fourier transformation to momentum space can be carried out. We 
do also a Fourier transformation of imaginary time to (fermionic) Matsubara 
frequencies Un = {2n+l)'7i/ (Ph) of the normal propagators (see AppendixEJ 

Gs{xt; x't') = Gs{x — x';t — t) 

= rn^l (0e— (-')e-(---')G,(fc,..„) . (1.17) 

with analogous definitions being valid for F{k, Un) and F\k, u;„). With these 
definitions, the propagators in Fourier space have dimensions of time. 

Introducing now = (^s to measure all excitation energies from the 
corresponding chemical potentials, the dynamical equations in Fourier space 
turn out to be 



{ihWn - ^k^)G^{k,UJn) + /^k{x)F\k,Un) = H , 
-ihujn - ^kl) F''{k,uJn) - Al{x)G^{k,uJn) = , 



{ihujn - ^fct) i^n) - Ak{x)Gi{-k, -UJ. 
[-ihjn - ^ki) Gi{-k, -Un) + Al{x)F{k, Un) = h , 



(1.18a) 
(1.18b) 
(1.18c) 
(1.18d) 



where we defined 



Ak(x) 



V 



d'yA{x,y)e'''- 



ik-{y—x) 



[1.19) 



We will consider in this work only translationally invariant systems. There- 
fore, the a;-dependence can be dropped: Ak{x) = A^. 
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Equations (I1.18|] form a system of four coupled, algebraic equations. To 
solve them, let us introduce the symmetric and antisymmetric (with respect 
to the exchange of t and | labels) quasiparticle spectra 

Es:=^^^ = e,-^^^e,-,^i,^ (1.20a) 
Ea ■■= ^^^y^ = -^^I^ ^ -51, . (1.20b) 

Here we have also defined the mean chemical potential /i = (/i| + /ij) /2 
and half the difference in chemical potentials = (/i| — /ij,) /2. With all 
these definitions, it is easy to see that Eqs. (I1.18cll.l8djl are just the same 
as Eqs. ()1.18b|1.18ap with the replacements Gi{k,u!n) —Gi{—k,—Un), 
A A*, F^{k,uJn) F{k,Un) and Eg -Eg. The solution is 



ihwn - Ea + E, 



s 



Gf(k, iUn) = h— Or. 

^ itfkUn-EA)'-El-\Ak\^ 

ifk^n + efc - (/i - (5/i) 

= h 2 (1.21a) 

{ihhJn + 5fi) - (efc - /i)^ - |Afc|2 



A 



k 



F(k, ujn) = —h- „ 

{inwn-EAf -El-\A, 

Afc 



-h 2 • (l-21b) 

(ifkJn + Sfl) - (Cfc - /U)2 - |Afc|2 



The solution for Gi{k,u!n) is identical to G'|(fc,c<j„) with the substitution 
Ea —Ea, while F'^{k,uJn) = [F{k,Un)]* . As expected, in the symmetric 
case (i. e., 6fi = 0) one can check that G'|(fc,u;„) = Gi{k,LJn). 

The quasiparticle excitation spectra of the system are just the poles of the 
propagators at zero temperature. These propagators can be found from those 
at finite temperature by analytic continuation of the imaginary frequencies 
onto the real axis iujn — > a; G R, and we have 

E^ = -6fi ± ^(efc-/i)V|Afc|2 = Ea± ^^i + |Afc|2 , (1.22) 

which gives the typical BCS dispersion relation E^ = ^/ E^ + | A^l^ in the 
symmetric case. In this case there is a minimum energy \Akp\ required to 
excite the system, from where the name '(energy) gap' derives. 

Substitution of the solution ()1.21bp into the definition of the gap func- 
tion ()1.19p yields the gap equation of a density-asymmetric system at finite 
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temperature: 



= j d^yV{y) [-F(Or+,yr)] e"='-^ 



'^k ,o , ,71,' A 



where we have identified the doubly-underlined term as the matrix element of 
the potential in fe-space, and used the techniques of Appendix 1X1 to evaluate 
the summation over Matsubara frequencies. 

1.2.3 An easy example: the symmetric case 

When both fermionic species have equal chemical potentials (which in the 
homogeneous system corresponds to equal densities), the normal propagators 
for up- and down-particles are the same and equal to 

with 



which are the typical results of BCS theory for quasiparticles with excitation 
spectra E^. can be interpreted as the probability amplitude of the paired 
state \k ],—k [) being unoccupied, and the probability amplitude of 
it being occupied, that is, there being simultaneously an up-fermion with 
momentum k and another down-fermion with momentum — fc, thus forming 
a Cooper pair. Then is the excitation energy of the (A^ + l)-particle 
system when one particle is added to the ground state of the A^-particle 
system, or the excitation energy of the (A^ — l)-particle system when one 
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particle is removed from the ground state of the iV-particle system |Sch88| . 
We show in Fig. 11.21 characteristic curves for m| and f |, together with Ek. 
From this plot the motivation for the name 'gap' is quite apparent: A is the 
energy distance between the ground state of the paired system and its lowest 
excitations: 



Ek\ 



A = A, 



In the general asymmetric case we have two different branches, one of which 
can be gapless; we will comment further on this point in chapter El 



Figure 1.2: BCS results 
for A/fx = 0.03: (Top) 
Quasiparticle distributions u| 
(black) and w| (green) and 
(bottom) quasiparticle energy 
Ek (black) compared to the 
ideal dispersion relation e^ — n 
(green) function of mo- 

mentum. The inset is a mag- 
nification of the region around 
the Fermi surface, which shows 
the existance of a minimum en- 
ergy value A required to excite 
the system. 




Momenlym k/k^^ 



The value of the gap in a symmetric system is also readily calculated 
starting from Eq. ()1.23j) . As now Ea = 0, we have E'^ = Ej^- 
ME^) - fpiE^) = -tanh(/?Efc/2). Therefore, 



E/^ , so that 



A./ 



7 — 7^ Vk'k — — tanh ( 3— 
(27r)3 " " 2Ek V 2 



(1.25) 



This equation is to be solved together with the number equation [compare 
with Eq. (ICTI) ]. 



d^A; 
(2^ 



— ^ ^ tanh ( 3— 
v/(efc-/i)2 + |A|2 V 2 



(1.26) 
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which accounts for the relationship between the chemical potential and the 
density of the system. However, in the weak-coupling limit in which we shall 
be working, one can safely take fi = ep. Then Eqs. ()1.25|) and ()1.26|) decouple 
and we only need to worry about the first one. 

There are two limiting cases of special interest: 
Zero temperature gap Agym 

When we set T = 1/(^^/5) = 0, the gap equation reads 

For a contact interaction of the form V{x — x') = g6{x — x'), we have 
Vfc'fc = 9 = constant, and the integral is ultraviolet divergent. In fact, this 
formulation of the gap equation has problems for potentials that do not have 
a well-defined transformation into Fourier space. It is better to resort to 
the T-matrix, that is the solution to the Lippmann-Schwinger equation and 
represents the repeated action (summation of ladder diagrams, see Fig. ll.Hjl 
of the potential V" in a propagating two-body system |Mes99| IGP89j 

j = y + y 1 v = V + V- ^ T, (1.27) 

E-H + ie E-Ho + ie ' ^ ' 

and is generally well-defined for any interaction. Then, the gap equation 
reads 

f d^k ^ f 1 1 \ , , , 



which is more convenient for numerial treatment (see Section E 

For dilute systems such as the atomic gases under current experimental 
research, one has kplal <^ 1. Therefore, as the integrand in the gap equation 
is sharply peaked around the Fermi momentum fci? in the weak-coupling 
regime (A//i <^ 1), we can substitute the momentum-dependent T-matrix by 
its value at low momenta, which is a constant proportional to the scattering 
lenght a. 

To := lim T^'k = • 1-29 

As To is a constant, it can be pulled out of the integral, and the gap function 
becomes independent of momentum. With this simplification, and making 
a change of variables to x = efc//i = /i^fc^/(2m/i) and d = A//i, the gap 
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Figure 1.3: Diagrammatic representation of the T-matrix: the sohd hnes 
represent free propagators, the dashed hnes stand for the bare potential V 
between the two particles, and the wiggle is the T-matrix. 

* 1 ^ -^-t 

! V + ! V ! V + ••• 
* 4 ► ^ 4— 




equation ()1.28|) can be rewritten as 



Tn 



47r2 



dkk'^ 



dx 



X 



X 



efc 



(1.30) 



where we have already integrated the angular variables and used the approx- 
imation fi ^ ef = h?kp/{2m) valid in the weak-coupling limit |FW7HlSch88[ 
IPB99| . From this result, we see that the gap in this limit will be a function 
only of the adimensional parameter kpa, irrespective of the specific form of 
the interparticle potential. Moreover, we shall have a non-vanishing gap only 
for an attractive interaction, a < 0, as in the case of the Cooper problem 
(Section OH). 

A quick, though approximate, evaluation of this integral can be done 
noting that, for (i -C 1, the integrand is sharply peaked around the Fermi 
surface (x ~ 1), which allows us to simplify the first term inside the square 
brackets taking a; = 1 in its numerator. The resulting integral is analytic, 
and we obtain 



isym 



(2^) 



A more careful evaluation of integral (ll.30jl gives the well-known result 
first obtained by Gor'kov and Melik-Barkhudarov |(;MR61I IPR99| 



A, 



sym 



/i 



1.083exp(-^) (a < 0) . (1.31) 

\2kFa J 



Note that the functional dependency on kpa is the same as in the previous 
expression and only the pre-exponential factor has changed. 
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Critical temperature Tc 

The critical temperature is defined as the lowest temperature where the gap 
vanishes. Assuming again a low density, so that the gap is independent of 
momentum and A^/ cancels with Afc in the numerator of the gap equation, 
we can find from (|1.23jl while setting /S. = Q m Eu'- 

^ = -9 f 7^7^tanhf/3,|~) ( 



This equation is solved in a way similar to that for the zero-temperature 
gap |GMB61] . The final result is 

where 7 = 0.5772 ■ ■ ■ is Euler's constant. Note that has the same func- 
tional dependence on kpa as the gap. Therefore, in the weak-coupling limit, 
the BCS theory leads to a simple proportionality relationship between the 
zero-temperature gap and the critical temperature for the superfluid transi- 
tion, independent of the physical system under study: 

^ = -- 0.567. (1.33) 

^sym 

This relationship — which is reasonably well fulfilled by many supercon- 
ducting metals and alloys |FW7H lKit96j . but not by high-temperature su- 
perconductors [CSTLOBt IBra98j — will allow us to calculate only the zero- 
temperature gaps. The transition temperatures into the superfluid state for 
dilute, atomic gases can be obtained then from Eq. (I1.33jl . 



1.3 Summary 

In this chapter we have introduced the pairing solution for the ground state 
of a system of interacting fermions. First, we have presented the concept 
of 'Cooper pairs'. Then, we have formulated the BCS theory as a solution 
in terms of Greens' functions at finite temperature of the interacting Hamil- 
tonian in a self-consistent Hartree-Fock-like approximation, in the general 
case of a two-component system with density asymmetry, i. e., where the 
two components may have different densities. 

Finally, we have obtained some analytic results for the gap at zero tem- 
perature and the critical temperature for the case where both components 
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have equal densities. In chapterOwe will study how the energy gap, chemical 
potential and other physical parameters of the sytem are affected when this 
last condition is not fulfilled. 

We remark at this point that an equivalent solution for the Green's func- 
tions for the full interacting Hamiltonian, can be found solving the corre- 
sponding Dyson's equation, which in (k, ci;)-space reads 

G{k,uj)-^ = G{k,uj)Q^ - S(fc,cj) . 

Here Gq^ = uj — + ii] sgn(c<j) is the non-interacting propagator and S is 
the self-energy, which has to be determined in order to solve the problem. 
A further analysis of this approach is beyond the scope of this introductory 
chapter and can be found in |SMPM05| . 
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Chapter 2 



Pairing in density-asymmetric 
fermionic systems 

Politiker, Ideologen, Theologen und Philosopher) versuchen immer 
und immer wieder, restlose Losungen zu bieten, fix und fertig gek- 
larte Probleme. Das ist ihre Pflicht — und es ist unsere, der Schrift- 
steller, - die wir wissen, dass wir nichts rest- und widerstandslos 
klaren konnen - in die Zwischenraume einzudringen. 

Heinrich Boll, Nobel lecture 

2.1 Introduction 

In Chapter [U we have developed a general formalism to study the BCS pairing 
mechanism in a system composed of two fermionic species with chemical 
potentials /i^ and /ij^, which are not required to be equal. This is a subject 
of interest in a number of fields of research, such as nuclear physics |MSn3bj . 
elementary particles |CNn4| and condensed matter physics |Yehn2| . We are 
mainly interested in the application to dilute, atomic systems where a great 
deal of work has been devoted since the pioneering contribution of Stoof and 
coworkers |^HSH96. H FS+QTj . who predicted the possibility to detect a phase 
transition to a superfluid state in a trapped system of ^Li atoms. 

After the achievement of Bose-Einstein condensation (BEC) in dilute 
gases of alkali atoms in 1995 |.TILA95[ IMIT95| . an important target in cool- 
ing atomic samples was to reach the degenerate regime in a gas of fermionic 
atoms. Indeed, the mechanism that allowed the production of BECs, known 
as evaporative cooling, relies intrinsically on the capacity of very cold atoms 
to quickly re-thermalize when some of them (those with higher energies) are 
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let escape from the trap, see Fig l2.ll At the very low temperatures of interest 
(of the order of fiK or below), the only remaining interaction in a dilute sys- 
tem are s-wave collisions. As these are forbidden between indistinguishable 
fermions by Pauli's exclusion principle |GP89[ lMes99j . this mechanism does 
not work to cool one-component Fermi gases {e. g., a gas with spin-polarized 
atoms). For such a system, the loss of very energetic atoms implies, of 
course, a decrease of the mean energy per particle, but in the absence of 
a (re)thermalization mechanism, this just means that the system is not in 
thermal equilibrium, but no redistribution of the remaining atoms in phase 
space occurs. 

Figure 2.1: Evaporative cooling procedure: in a trapped system of bosons 
at temperature T, the most energetic ones are let escape, thus reducing the 
mean energy per particle of the remaining ensemble. After rethermalization, 
the system is at a lower temperature T' < T. 




T <E'> < <E> T'<T 

This limitation can be overcome if the trapped system is not composed 
just by one kind of fermions. For example, one can trap atoms in two (or 
more) hyperfine states. In this case, there is no problem for the existence of 
s-wave collisions between atoms belonging to different states and this simul- 
taneous cooling scheme works for such a multi-component fermionic system 
in the same way as for a bosonic system, as was first shown by DeMarco 
and Jin at .JILA |Rice99j . Another possibility is to trap a mixture of bosons 
and fermions. In this case, the bosons are cooled in the usual way, while the 
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fermions cool down by thermal contact with the bosons, as s-wave collisions 
between them are not forbidden. This mechanism is known as sympathetic 
cooling and it was first realized by the group of Randy Hulet at Rice Uni- 
versity, who reported experiments where a gas of ®Li atoms (with fermionic 
character) was driven to quantum degeneracy by contact with a gas of ^Li 
(bosons) |Riceni| . 

Since these pioneering works, many groups have produced degenerate 
Fermi gases around the world [JILA, MIT, LENS, Innsbruck, Paris,...]. One 
of the main goals of this research has been to create (and detect!) a super- 
fiuid made up of atoms, in a sense analog to superfiuid helium but with the 
advantage (in principle) of a much weaker interaction due to the low density 
of these systems (typically, in the absence of resonant phenomena such as 
Feshbach resonances, p\a\^ ^ 1, while in helium pa^ ~ 1), which allows for 
an easier understanding of the physics behind the experiments. In particular, 
it is possible to use the mean-field BCS theory of pairing in the weak-coupling 
approximation |Fj1HSH96[ IhFS+97| . 

In this Chapter, we will study which are the possibilities of having pairing 
in a two-component fermionic system such as that in |Rice99j . In particular, 
we will focus on the relevant paper played by the difference in densities of 
the two components. As we will show, this difference reduces dramatically 
the size of the energy gap (and, therefore, the expected critical temperature 
for appearence of superfiuidity) when considering pairing between atoms be- 
longing to different species. However, pairing between atoms belonging to 
the same species can be enhanced by properly adjusting the density of the 
other species. 

2.2 Pairing in s-wave 

Let us assume a fermionic system composed of two distinct species (which we 
label t and |) of equal mass m and with densities and pj, or equivalently 
total density p = Pt + Pi ^^'^ density asymmetry a = {pi — Pi)/ {p] + Pi)- 
For example, species t and | can denote different hyperfine levels of atomic 
gases, or 'neutron' and 'proton' in the context of nuclear physics (where one 
should take into account both the spin and the isospin as internal degrees of 
freedom of the 'nucleon'). We also introduce the Fermi momenta ks = k^^ = 
(Gvr^Ps)^/^, and chemical potentials ps = h?k'l/{2m) {s =t, |), together with 
p = (pt + /ii)/2 and 6p = (pt - p^)/2. 

For simplicity, we will for the moment consider an idealized system where 
t and I particles are interacting via a potential V with s-wave scattering 
length a, while direct interactions between like particles j-j and |-| are 
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Figure 2.2: The two possible lowest-order pairing interactions: (a) Direct s- 
wave interaction between different species, (b) Polarization-induced p-wave 
interaction between like species. Vq and Tq are the s-wave (L = 0) bare 
potential and T-matrix between species | and I, respectively. 




(a) (b) 

absent [cf. Eq. ()1.4|1 and Fig. I2.2f a)]. We are interested in the situation at 
very low density, ks\a\ <C 1. We have seen in Chapter [T] that, in this limit, 
the pairing properties of the symmetric system are completely determined by 
the scattering length or, equivalently, the low-momentum s-wave T-matrix, 
To = Anfi^a/m. We will now study the general case where the two species can 
have different densities, and analyze the corresponding pairing gaps generated 
by this interaction V. 

As we did for the symmetric system in Sect. I1.2.3| we must solve the gap 
equation (ll.23jl . 

^^'-/|^^-^[^-(^^')-^-(^^")]' 

E^ = Ea± y^^| + |AfcP , 

together with the equations that fix the chemical potential of each species (or, 
equivalently, fi and 6fi) from the values of their densities. These equations 
can be found again starting from the normal propagators and summing over 
Matsubara frequencies (see Appendix : 

= <fiE^) + vlfiK) = ul {]{El) - fiE,)] + f{E^) , (2.1a) 
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These distributions are shown in Fig. 12.31 for = 0.1 and different values 
of the gap A to get some physical insight. The sharp Fermi surfaces of the 
non-interacting system smear out for increasing A reflecting the correlations 
introduced by the pairing interaction. However, this 'smearing' hardly pen- 
etrates into the region between /ij and /i^ as Pauli's principle forbids new t 
particles to enter there. Thus, the newly formed j-j pairs need to climb to 
states efc > This has a cost in kinetic energy, which is compensated by 
the pairing energy gained. Clearly this mechanism cannot sustain arbitrarily 
large asymmetries, as the kinetic energy investment grows rapidly with the 
width ~ 25 e of the forbidden region (see below). In fact, one can readily find 
the total density p and density difference 5p to be 



P 



k ^ 

5P = PT - Pi = E [f^^^t) + fF{E^) - 1] . (2.2b) 



At zero temperature, friE) = Q{—E) and, therefore, 

ME^) + ME^) - 1 = 0(-E+) , 

where we used the fact that 6(x) = 1 — 6(— x). From the second equation 
we see that the unpaired particles are to be found in the energy interval 
[fi — 6e, p + 6e\, with 6e = a/ — A^ (^ Sp for A < p), which does not 



Figure 2.3: BCS momentum dis- 
tributions of the species | (black 
lines) and | (green lines) in asym- 
metric matter. The continuous 
lines stand for the normal mat- 
ter (A = 0) distribution, the 
dashed lines for superfluid matter 
with A = 0.03/i, while the dotted 
lines show the distributions cor- 
responding to the more coupled 
case A = 0.07/i. [see Eqs. (l2lTl ]. 
The shaded area is the 'forbidden 
region' [p — 6e,fi + 6e] for the case 
A = 0.07/i. 
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contribute to the pairing interaction, as indicated by the shaded area in 
Fig. \2.'M This leads to a rapid decrease of the resulting gap when increasing 
the asymmetry.* 

In the weak-coupling case, A ^ (5e <^ /i, which is adequate in the low- 
density limit, the momentum distributions of the two species are anyhow 
very sharp and one obtains from Eqs. (I2.2jl 



1 + f 



(27r)3 

1 /2m/i\^/' ^ ^ n+s. 



1-^^ x-1 r , ^ x-i 



dx \fx — -j^^^=^^^= — / dx \fx 







27r2 V ^2 / 3 



and 



bp 



47r2 



^ (^) ^x^Q (5x - v/(x - 1)2 + ci2 



27r2 V 

where we used the shortening 5x = Se/n and d = A//i, and the fact that 
6fj, <ti p. Therefore, 

a = ^ « -- « 1 , 2.3 
p 2 p 

i. e., the width of the forbidden region 6e is directly proportional to the 
density asymmetry a. Analyzing in a similar way the gap equation, one 
can obtain these relationships involving the parameter Se jSHSHQBj ISAL971 

EEnni, 

Sn + Se = const. = A^y^ ^ = A^y^^ - 2A^y^5e , 



*One could imagine also a promotion of all the t particles responsible for the block- 
ing effect to higher kinetic energy and having a continuous distribution of pairs at low 
momenta, but this would imply a larger kinetic energy investment. 
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where Asym is the gap in symmetric matter of the same density p [Eq. JUSH)]- 
Altogether, we can write the gap as a function of the asymmetry thus: 

It vanishes at a^Sf = 3Asyin/(4/i), which is an exponentially small number 
in the weak-coupling limit kp\a\ <^ 1. Therefore, for very small asymmetries 
already, pairing generated by the direct interaction between different species 
becomes impossible according to the BCS theory. However, this limitation 
can be partially overcome, as we will see in chapter |H1 



2.3 Numerical results 
2.3.1 The algorithm 

Let us now present our numerical procedure to solve the gap equation and 
check the results with the previous analytical calculations. We start rescaling 
all energies in units of the mean chemical potential, 

x = eu/p d = A//i 5x = 5e/p 

= Ejfi E^ = -^± E, = E^/fi , 

Assuming a constant value for Vfc/fc = g, so that the gap function will not 
depend on momentum, we rewrite the gap equation as 



1 = 9 



y/JI dx [fp (px^) - fp {fix')] . 

Jo 



To avoid a possible divergence in the gap equation due to the use of a contact 
potential, we resort to the T-matrix as in Sect. 031 (see also |MPS98j ). We 
get for the general low-density, asymmetric case 



where we used To = inh'^a/m for a low-density system. Regarding the 
densities, they are given by integrating Eqs. (I2.H1 . We present below our 
numerical results for the solution of equations ^2.1^ and ()2.5|1 obtained in the 
following way: 
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1. Introduce an 'effective Fermi momentum' fc^ through the chemical po- 
tential, /X = fi^/c^/(2m). 

2. Define the functions 

^ " ~ / da;yx{l - frifJ'X'^) + ul [/^(/ia;^) - } 

with ul = [1 + (x — l)/ Ex]/2. These are functions of the gap d (through 
Ej.) and the chemical potentials /i| and /i^ (through /x, /c^ and x^). 
They must be zero at the solution of Eqs. (I2.H1 and ()2.5jl . 

3. Assume some initial values d = d^^^\ fi-^ = /ij^^ As we work at fixed 
total density p and density asymmetry 6p = pa, we have used these 
guesses: 

d^'' = ^ [cf. Eq.JZl], pf^=e,{l±af\ 
p 

with ep = ^2(37rV)^^V(2r«). 



4. Insert these values into / and find (e. 5^., by a Newton- Raphson method) 
a zero of / as a function solely of d. Call this zero d^^\ 

5. Insert d^^^ and p^^^ into g, and find a zero of 5^ as a function of solely 
p^. Call this zero p^^\ 



6. Insert d''^\ p^^^ and into h, and find a zero of /i as a function of 
solely pi- Call this zero p^^\ 

7. If the new values for d and pn make |/|, l^fl and \h\ smaller than 
a certain tolerance e (in our calculations, typically e ~ 10~^), take 
them as the solution. Otherwise, go back to step 4 and try again until 
convergence is reached. 
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2.3.2 Symmetric case 

First of all, in order to test our code, we analyze the dependence of the gap 
and the chemical potential for the symmetric case, k-^ = = kp, at low 
temperatures. We fix T = O.lTc, with Tc the BCS weak-coupling prediction 
for the critical temperature as given by Eq. (jl.HHj) . The results that we obtain 
are plotted in Fig. 12.41 In the left panel we show with circles the calculated 
variation of the energy gap A with the coupling fci?|a|. For comparison, the 
solid curve is the BCS prediction at low-density, Eq. (|1.3H1 . We see that 
the agreement is very good for all kpla] < 1, and only deviates slightly for 
/ci?|a| > 1.3. The 'effective' Fermi momentum fc^ = ^/Q/mJI/h (normalized 
to the non-interacting Fermi momentum kp = y/2meF/h) is plotted in the 
right panel also as a function of the coupling. In this case, we find A;^ ~ kp 
for kp\a\ < 1, with a smooth decrease that becomes more pronounced above 
this value of the coupling. This is a signal of the increasing modification of 
the system with respect to the non-interacting one for such large couplings. 
In fact, for these values of kp\a\ the mean-field BCS theory is not adequate, 
and more sophisticated methods need to be used |HKCWni| 10002^ IBYn4[ 
IFSn4[ lPPSn4[ ICSTLfl5j . Therefore, we shall not work in the following in this 
strong-coupling regime. 

Next we analyze the temperature dependence of the same quantities 
above. To satisfy the weak-coupling condition just discussed, we fix the 
density to p = 2 x 10"^^ cm"^, so that kp\a\ = 0.445524797 < 1. The results 

Figure 2.4: Solution of the BCS equations for a temperature T = O.lTc 
and varying coupling: (Left) Numerical solution for the gap A as a function 
of the coupling and prediction of BCS at low density, Eq. (ll.Hljl . (Right) 
Effective Fermi momentum k^ as a function of the coupling. 
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are shown in Fig. l2.5[ In the left panel, the energy gap [divided by the zero- 
temperature value (IHSH)] is plotted as a function of the reduced temperature 
T/Tc, with Tc given again by Eq. ()1.33|) . Our numerical results are shown 
with empty circles, while the two lines are the analytical curves of the BCS 
model (see Ref. |FW71j . p. 449) at low and 'high' temperatures: 



Mow T — ^sym 



1 _ Javr^e-^--/'^-^ 

^sym 



^high T 



I. T (^ ^ 



1/2 



T/T, < 1 
(T, - T)/T, < 1 



Both predictions agree well with the data in the corresponding ranges of 
validity. More specifically, the low-temperature curve is very close to the 
calculated values for T/Tc < 0.5, while the 'high'-temperature one is valid 
for T/T, > 0.85. 

Regarding the chemical potential (shown by the circles and the solid 
curve, which is a guide to the eye), it is always very similar to the non- 
interacting one, i. e., Sp- However, the weak attraction between the fermions 
results in a slight reduction of /i, that diminishes when approaching the 
critical temperature, where the gap is expected to vanish (but the system is 
still interacting, hence the difference between /x and Sp). The dependence of 
jjL on temperature can be well fitted by 

/i = /i(T = 0) [1-7(T/T,)^] , 

Figure 2.5: Solution of the BCS equations at finite temperature for a fixed 
density p = 2 x 10~^^ cm~^: energy gap (left) and chemical potential (right) 
vs. reduced temperature T/Tg. See the text for details on the notation. 
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as shown by the dashed curve in the figure. Here /i(T = 0) is the chemical 
potential at zero temperature, and 7 an adjustable parameter, for which we 
found the best value 7 = 4.28 x 10"'* (xlt = 3-06 x 10"^°). In an equivalent 
plot, the inset to the right figure shows the variation of l — k^/kp as a function 
of temperature. 

2.3.3 Asymmetric case 

After this short overview of the symmetric case, we will study the dependence 
of the gap on the density asymmetry a = 5p/ p in order to check numerically 
Eq. (I2.4jl . We consider T = and again the density p = 2 x 10^^ cm~^ 
(fc^lal ~ 0.4) as we have seen in the previous paragraph that it is well into 
the BCS validity regime. Our results are summarized in Fig. 12.61 On the left 
panel, we show how the energy gap A varies with increasing asymmetry a. 
The circles are the calculated data, while the line corresponds to Eq. (j2.4j) . 
with the Asym given by Eq. the dashed vertical line indicates the 

expected asymmetry a^^^ = 2.4% where the gap will vanish. The agree- 
ment between the calculation and the expected behavior is excellent for all 
asymmetries, even very close to this limiting asymmetry. 

On the right panel, we show, for comparison with the symmetric case, 
the dependence of the mean chemical potential p together with the chemical 
potentials p^ and 6p. We see that p depends only very weakly on a, as 
compared to what happens with p^ and pi- These start being equal in the 

Figure 2.6: Solution of the BCS equations for a fixed density p = 2 x 10~^^ 
cm~'^ at zero temperature as a function of the density asymmetry a. (Left) 
Energy gap; (Right) Chemical potentials: total p (circles), p^ (triangles up), 
Pl (triangles down) and difference 6p (filled squares). 
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symmetric case, and separate in an approximately linear way for increasing 
a until a ~ 0.01, where they 'saturate' to the values that they will have for 
a = a^^^. As one would expect, an analogous behavior is followed by 

As a conclusion, we have checked numerically the expected behavior of 
the gap for the s-wave solution to the coupled equations determining the 
gap and the densities. In particular, we have seen that the mean chemical 
potential is essentially equal to the non-interacting value ep and that the gap 
vanishes for very small asymmetries. At finite temperatures, as one would 
expect, we have checked also that the gap decreases, and so the maximum 
asymmetry that would allow pairing is even smaller. 



2.4 Pairing in p-wave 

Therefore, for larger asymmetries only pairing between identical fermions can 
take place. In this case, Pauli's exclusion principle demands the paired state 
to be antisymmetric under exchange of the particles, which for spin-polarized 
atoms means that they must be in a state of odd angular momentum L. The 
leading contribution to the interaction will be the p-wave one. 

We discuss in the following the pairing of j-particles mediated by the 
polarization interaction due to particles pertaining to species | as shown in 
Fig. I2.2r b). We will check below that this contribution is dominant over 
the direct p-wave t-T interaction at low density. Quantitatively the relevant 
interaction kernel reads |(;MH6H IHRSVOOj ISPU01I IK(]88[ IBKK96| 



{k'\T^\k) 



T 



(2.6) 



The factor 1/2 corrects for the fact that conventionally the Lindhard function 
(pertaining to species |), 



mk\ 



1 1 

2 ^ 



X 



In 



1 + x 



1 



X 



X 



(2.7) 



contains a factor two for the spin orientations |FW71| . which is not present in 
our case as we are dealing with fermions without internal degrees of freedom. 

One can see that IIj^ < 0. Therefore, the effective interaction is always 
attractive, irrespective of the sign of a |FL68] . This is due to the absence of 
exchange diagrams, which makes IIj^ attractive in contrast to the case of one 
species, where the polarization effects reduce the s-wave BCS gap by a factor 
(^4e)i/3 ~ 2.2 even in the low-density limit, where one would naively expect 
them to be negligible. |(;MR61 I IHPSVOOL ISFRlVl] . One can therefore assert 
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that, even in the limit kp ^ BCS is not such a good aproximation, as it 
disregards these higher-order interactions. 

Now we project out the p-wave («. e., L = 1) contribution of the inter- 
action with the Legendre polynomial Pi{z) = z, with z = k' ■ k the cosine 
of the angle between the in- and out-going relative momenta, which in the 
range of the weak-coupling approximation we can take on the Fermi sur- 
face |KC89l lEMRKOOl 



.{L=l) 



1 

' ^ 
2 2 



j ^ dzzli^{^2{l-z)k^). 



The integral of the Lindhard function can be evaluated making use of the 
result |SPRnH 



dxx Yiix) = — 
^ ' 12 



2\n\l- Z^\ + {5-3Z'^)ZHn 



1 - Z 



+ 2Z^ + 6Z^ 



resulting in 



a2A;^21n2-l /A; 



m 



-9 



k, 



where we introduced 

9{y) - 



6(2 In 2 - 



(4-102/2)ln|l-y2 



l + y 



■ (2.8) 



This function is normalized in symmetric matter, g{l) = 1, and it is plotted 
in Fig. I2.7r a.). It is customary to replace the numerical factor (2 In 2 — l)/5 
in Eq. (I2.4|l by its approximate value 1/13. 

We use now the general result for the (angle-averaged) L-wave pairing 
gap iFLBlfMTml lKLBBl. 



AL(fct) CL/i^exp 



2n^h^ 



A;^Tl(%,/c|)_ 



(2.9) 



with Tl the L-wave projection of the relevant interaction and cl a constant 
of order unity. Considering that for a pure polarization interaction F^, the 
leading order in density is = F^,, we get 



Ai(/c|) = Cl— exp 
2m 



13{7i/2f 



a?k^k^g{k^/k^ 



(2.10) 
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Taking into account the dependence of this expression on the two Fermi 
momenta = [37r^p(l ± a)]^^'^, the final result for the variation of the 
pairing gap with asymmetry a and total density p can be cast in the form 

1^ = {l + ar/'exp[u{p)h{a)] (2.11) 



with 



1 

^^"^ " ^ ~ (l-a2)i/3^[((l + «)/(l-a))i/3] (2-^2) 



and the density parameter 

2 



The function h{a) is displayed in Fig. I2.7r b). where a = 1 corresponds to 
pure t-matter. One notes a maximum at Oopt ~ 0.478 with an expansion 
h{a) ^ 0.465 - 1.343(a - 0.478)^. 

At the optimal asymmetry, the gap is enhanced with respect to that of 
the symmetric case {h{0) = 0] by a factor e"-^^^" a; 10°'^°". In the low-density 
limit u ^ oo this represents an enormous amplification at finite asymmetry. 
Around this peak, the variation of the gap with asymmetry is well described 
by a Gaussian with width a = 1/(1. 64y^). As an illustration of this effect, 
Fig. l2.7r c) shows the ratio (|2.11|] for a value of the density parameter u = 100 
(corresponding to kpla] ^ 0.57). As expected from the previous discussion, a 
peak of the order of 10^° is observed, that becomes rapidly more pronounced 
and narrower with decreasing density p (increasing m), although of course 



Figure 2.7: (a,b) The functions g and h, appearing in Eqs. (I2.8jl and (j2.12|l . 
respectively, (c) The variation of the gap with asymmetry, Eq. (I2.1H1 . for a 
density parameter u = 100. 
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Figure 2.8: Third-order diagrams in a one-component system: (a) Direct 
j9-wave interaction, (b) Polarization contribution. 
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at the same time the absolute magnitude of the gap decreases strongly with 
decreasing density: Ai(p, 0) oc exp(— m) [Eq. ()2.10|) ]. 

Let us now briefly discuss higher-order efl'ects, namely contributions of 
order {kpa)^ the denominator of the exponent of Eq. (I2.1()jl . There are two 
principal sources of such effects, which are shown diagramatically in Fig. 12.81 
The first one. Fig. \2.H( a). is the direct p-wave interaction |YM99| between 
like species that we have neglected before. Parametrizing the low-density 
p-wave T-matrix in the standard form Ti{k,k') ATiafkh' /m, where af is 
the p-wave scattering volume, leads to a BCS gap |YM99| 



Ai 
/it 



=8/3 



exp 



TC 



2{k^a 



\3 



(2.13) 



where we have explicitly determined the prefactor of the exponential. This 
gap is much smaller than the |-mediated one for {|a|, |ai|} ^ k^^, which 
justifies the fact that we neglect this effect in the previous calculation. How- 
ever, around ap-wave Feshbach resonance, where |ai| ^ |a|, this effect might 
be important, even though this kind of pairing has not been observed in the 
experiments that explored the p-wave resonance of in the hyperfine state 
\F = 9 /2, mp = -12) |.TTT:An3c| or ^Li in \F = 1/2) |ENSr)4h| (see also 
[BohOO] for theoretical calculations of the position of the resonances, and 
[ .StanOOf Ilnnsn3| for experiments on bosonic cesium). 

The second type of third-order contributions are polarization effects in- 
volving the s-wave scattering length a. In contrast to the case of a one- 
component system, where there are several relevant diagrams, in the present 
two-component system only one diagram exists, Fig. I2.8r b). Unfortunately 
it can only be computed numerically, which we have not attempted. 
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Finally, to fourth order, there is a large number of diagrams contributing 
to the interaction kernel. Moreover, at this order it is also necessary to 
take into account retardation effects, i. e., the energy dependencies of gap 
equation, interaction kernel, and self-energy need to be considered [AG92]. 
All this can only be done numerically and was performed in Ref. (EMRKOnj 
for the case of a one-component system. 

In any case, the existence of higher-order corrections will not alter the 
main conclusions drawn so far, namely the presence of a strongly peaked 
Gaussian variation of the gap with asymmetry. They may, however, shift this 
peak to a different density-dependent location, and also modify the absolute 
size of the gap. Furthermore, we must note again that the perturbative 
approach that we have followed is clearly limited to the low-density range 
kp\a\ < 1. For larger couplings, different theoretical methods are required, 
which is an interesting field of current investigation |BB73[ IAB87[ IHKCWflT| 
Kxm iRYfMl i PPSlMl iFDSlMl K]STLn5| . 

2.5 Summary 

In this chapter, we have studied the possibility of pairing in a system com- 
posed of two distinct fermionic species, assuming that the direct p-wave 
interaction between like species can be neglected. First, we have studied 
analytically the gap produced by a direct, attractive, s-wave interaction be- 
tween different species. We have shown that it produces a gap ~ exp [7r/2/c j?a] 
[Eq. JUSlI)] only for very small asymmetries between the densities of the two 
species, a < Asym//i, Eq. ()2.4p . We have numerically checked this behav- 
ior, and we have seen that, in fact, the maximum asymmetry decreases with 
temperature. 

However, a p-wave attraction between two fermions of the same species 
produced by a polarization of the medium of the other species, can give rise 
to pairing over the whole range of asymmetry. In practice a sharp Gaus- 
sian maximum around a ~ 0.478 {p]/pi ~ 2.83, k^/ki ~ 1.41) appears. 
Unfortunately, the absolute magnitude of this gap is rather small, as it de- 
pends strongly on the density of the system, ~ exp[— 13(7r/2A;i;'a)^], which 
we assumed low in order to be in the range of validity of the BCS mean-field 
approach. 

We argue that higher-order corrections may modify quantitatively but 
not qualitatively these general features. 

The experimental observation of both types of pairing in dilute, atomic 
vapors is expected to be difficult. For the case of p-wave pairing due to 
the extremely small size of the expected gap, which translates into a very 
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low transition temperature [remember Eq. (jUHBl)], which is difficult to reach 
experimentally (see Sect. 12. ill . For the s-wave case, the nearly perfect sym- 
metry that is required poses a difficult experimental challenge. However, the 
present possibility to transform a two-component Fermi gas into a molecular 
gas by crossing a Feshbach resonance |,nLAfl4b| IInnsn4j opens the door to 
create a perfectly symmetric system by removing the atoms not converted 
into molecules and, then, going back across the resonance, forming again a 
two-component fermionic system, now with well-balanced populations. Up 
to now, the experiments have focused on the strongly-interacting regime 
kplO'l ^ 1 where superfluidity might arise at higher temperatures. Once 
this has been succesfully achieved |MITn5| . one could imagine to go further 
to the weak-coupling regime and transfer atoms in a controled way from one 
hyperflne state to the other by RF pulses, which would allow a precise exper- 
imental check of Eq. ()2.4p in case low enough temperatures were achieved. 
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Chapter 3 



Pairing with broken space 
symmetries: LOFF vs. DFS 




Bill Waterson, Calvin and Hobbes 

3.1 Introduction 

In chapter n] we introduced the BCS theory for pairing in a fermionic system. 
We showed that this theory predicts a phase transition into a superfluid state 
for temperatures of the order of Asym/fc^, where in the weak-coupling regime 
krlal < 1 [see Eq. ^TTH ] 

Stoof and coworkers were the first to point out that this result, when applied 
to dilute systems of ^Li, predicts relatively large gaps for moderate densi- 
ties (e. (/., p ~ 10^^ cm~^) due to the large (negative) value of the s-wave 
scattering length of this species, a = —21600^ (as ~ 0.53 A is Bohr's ra- 
dius). These large gaps would translate in transition temperatures of the 
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order of nanokelvin and, therefore, possibly within experimental reach in a 
short time after their proposal in 1996 [SHSH96] . Regarding the problem on 
the density asymmetry between the two species to pair, they also pointed out 
'that the most favorable condition for the formation of Cooper pairs is that 
both densities are equal'. Then, they solved the gap equation in this most 
favorable situation, for a trapped system of ^Li atoms using local density 
approximation (lda) [HFS+97J. 

In chapter Owe have carefully analyzed why a two-component symmetric 
system will always have s-wave gaps larger than asymmetric systems. Be- 
sides, we have also shown that these s-wave gaps are non-vanishing, within 
standard BCS theory, for quite a narrow window of density asymmetries a. 
For a > a^Sf = 3Asym/(4/i) [cf. Eq. (12. 4|] ]. only p-wave pairing between 
atoms in the same hyperfine state seems possible. However, we have seen 
that the size of this gap is rather small, and in fact it is difficult to expect 
this phase to be detected experimentally. Now we will explore another way 
of overcoming the density-asymmetry problem. To this end, we solve the gap 
equation on a wider space of functions, allowing for more complex structures 
than the typical rotationally-symmetric order parameter found in chapters [T] 
and El 

In references |MSn2[ IMSfl8b| it was shown that for nuclear matter at 
saturation density — which is a strongly coupled system (Ao//i ~ 0.3) — a 
superconducting state featuring elliptically deformed Fermi surfaces (dps) 
was preferable to the spherically symmetric BCS state. Following these ideas, 
we propose here to explore the weak-coupling regime of this theory in con- 
nection with current experiments with ultracold, dilute systems. 

It has been known for a long time that the homogeneous BCS phase can 
evolve into the Larkin-Ovchinnikov-Fulde-Ferrell (loff) phase [L064.. FF64j . 
which can sustain asymmetries a > ttmSf by allowing the Cooper pairs to 
carry a non-vanishing center-of-mass momentum. Note that, even though 
both the LOFF and the DFS phases break the global space symmetries, they 
do it in fundamentally different ways: the LOFF phase breaks both rota- 
tional and translational symmetries due to the finite momentum of the pairs' 
condensate. On the contrary, the DFS phase breaks only the rotational sym- 
metry [from 0(3) down to 0(2)]. We remark also that, after so many years 
since its proposal, only very recently some experimental evidences of the 
detection of the LOFF phase have been reported in the heavy-fermion com- 
pound CeCoIns (see |RFM+03t lRMC+n3[ IWKT+fMl IKSK+OBI IM AT+n5j V In 
this sense, it is interesting to note that atomic systems offer a novel set- 
ting for studying the LOFF phase under conditions that are more favorable 
than those in solids (absence of lattice deffects, access to the momentum 
distribution in the system through time-of-fiight experiments) as explored 
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in |(]omnH l(M02[ IMC ^031 IMMT05L IYRn05| . These systems also offer the 
possibility of novel realizations of the LOFF phase which for example in- 
voke p-wave anisotropic interactions [Comfllj . There has been also much 
interest in the LOFF phase in connection with hadronic systems under ex- 
treme conditions where the interactions are mediated by the strong force, 
see I ARRni L ISeHnTI IRRn2[ f(:;Nn4| . In this context no experimental detection 
of this phase has been reported yet. 

We note that another possible configuration for a density-asymmetric 
system is the phase separation of the superconducting and normal phases 
in real space, such that the superconducting phase contains particles with 
the same chemical potentials {i. e., it is symmetric), while the normal phase 
remains asymmetric, see [Bed 021 lB( ]Rn3| R ]ain4| . The description of such an 
heterogeneous phase requires knowledge of the poorly known surface tension 
between superconducting and normal phases, and will not be attempted here. 

We shall compare below the realizations of the DFS and LOFF phases in an 
ultracold gas of ^Li atoms where the hyperfine states |t) = \F = 3/2, mp = 
3/2) and ||) = |3/2,l/2) are simultaneously trapped. Here, F denotes the 
total angular momentum of the atoms in units of h, while mp is its projection 
onto the z axis of some reference frame. Fermionic systems where two hyper- 
fine levels are populated have been created and studied experimentally with 
^Li and ^"K atoms (see [RI^^ IRENS02L IDuke02[ IrII^ IMTT03L lENS03j . 
to quote a few examples). The mixture of different hyperfine components al- 
lows one to overcome the problem of cooling fermions set by Pauli's exclusion 
principle, as indicated in Sect. 12. ll 

These systems are characterized by a hierarchy of length scales. The 
largest scale is usually set by the harmonic trapping potential. As it is much 
larger than any other scale in the system, we will neglect the finite-size effects 
for the moment, and assume an homogeneous system for our analysis. A step 
further would be to perform an LDA calculation, in a way similar to that used 
in |HFS"'"97j for the standard BCS treatment or |CKMHfl2] for pairing in a 
resonantly interacting Fermi gas. 

The typical range of the interatomic, van der Waals forces is R < 10~^ cm 
while the de Broglie wavenumber of particles at the top of the Fermi sea is 
kp ~ 10^ — 10'' cm~^. Therefore kpR <^ 1, and the interaction can be 
approximated by a zero-range force characterized by the s-wave scattering 
length a. For the particular case of collisions between ^Li atoms in the above- 
mentioned states, a = —21600^, and we obtain kp\a\ ~ 0.04. Therefore the 
system is in the weak-coupling regime, since u{kp)g = 2kp\a\/'K <^ 1, where 
^{kp) = mkp/ {2TT^h?) is the density of states at the Fermi surface of the 
non-interacting system, and g = iiih'^a/m is a measure of the strength of the 
contact interaction. For larger values kp\a\ > 1, the bound states need to 
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be incorporated in the theory on the same footing that the pair correlations, 
and the formal treatment becomes more delicate. 

To study the LOFF phase, we will assume that its order parameter has 
a simple plane-wave form. Even though more complex structures for the 
LOFF order parameter can be studied, we believe this would not change 
qualitatively our results. Furthermore, we will show that, allowing the Fermi 
surfaces of the species to deform into ellipsoids, the range of asymmetries 
over which pairing is possible is enlarged with respect to the predictions of 
the standard BCS or LOFF theories. Also, the gap in asymmetric systems 
where pairing was still possible within those frameworks, will be shown to be 
larger in the DFS phase. 

Finally, at the end of the chapter, we shall describe an experimental sig- 
nature of the DFS phase that can be established in time-of-flight experiments 
and that would allow one to distinguish the DFS phase from the competing 
phases. 



3.2 Breaking the symmetry: LOFF and DFS 



3.2.1 Description of the LOFF state 

While the BCS ground state assumes that the fermions bound in a Cooper 
pair have equal and opposite momenta (and spins), for fermionic systems 
with unequal numbers of spin up and down particles this is not always true. 
In this situation, Larkin and Ovchinnikov |L()64| and independently Fulde 
and Ferrell |FF64| noted that the pairing is possible amongst pairs which 
have finite total momentum with respect to some fixed reference frame. The 
finite momentum P changes the quasiparticle spectrum of the paired state. 
To see this, we can write down the normal propagator in that reference frame: 



G'f^(fc,P;a;. 



{ p 



n -1 



(3.2) 



Now, using Eq. ()1.2n|l . the symmetric and anti-symmetric parts of the quasi- 
particle spectrum read 



Ea = ti- 



8m 2 
P k _ /it 

2m 2 



(3.3a) 
(3.3b) 



Fortunately, the results of the previous chapters remain valid with the above 
redefinitions of Es and Ea- Note that the quantities of interest, in particular 
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the gap, now depend parametrically on the total momentum. Interestingly, 
Ea in does not vanish in the limit of equal number of spin-up and down 
particles {i. e., when = /ij^). In this case, the LOFF state (a condensate of 
pairs all with momentum P) lowers the energy of the system with respect 
to the normal (unpaired) state. Nevertheless, it is not the real ground state 
of the symmetric system, as it is unstable with respect to the ordinary BCS 
ground state. In fact, it is well know that, for the symmetric system, the 
most favorable configuration for pairing is for P = 0. 

3.2.2 Description of the DFS state 

We now turn to the deformations of the Fermi surfaces. The two Fermi 
surfaces for spin-up and -down particles are defined in momentum space for 
the non-interacting system by the fact that the energy of a quasiparticle 
vanishes on them: 

When the states are filled isotropically within a sphere, the chemical poten- 
tials are related to the Fermi momenta ks (s=t, i) as /i^ = li?k'l/{2m) (for 
simplicity we assume here that the temperature is zero). To describe the 
deformations of the Fermi surfaces from their spherical shape we expand the 
quasiparticle spectra in spherical harmonics = ^i^s^ Pi{x)i where x is 
the cosine of the angle formed by the quasiparticle momentum and a fixed 
symmetry-breaking axis; Piix) are the Legendre polynomials. The / = 1 
term breaks translational symmetry by shifting the Fermi surfaces without 
deforming them; this term corresponds to the LOFF phase and is already 
included by using ep/2±fc- Truncating the expansion at the second order 
(/ = 2), we rewrite the spectrum in the form | |MSn2[ lMSfl8bj 

where the parameters ris = ^i^V/^s describe a quadrupolar deformation of the 
Fermi surfaces. It is convenient to work with the symmetrized S = {ri'^+rii)/2 
and anti-symmetrized 6e = (r/^ — r]i)/2 combinations of rj^^^. For simplicity, 
below we shall assume S = 0, and consider only two limiting cases: 5e 7^ 
and P = (the DFS phase) and 6e = and P 7^ (the plane-wave LOFF 
phase) [see Table EH] for a summary of the nomenclature used]. Thus, in the 
DFS phase we have 

^fcT/i = efc -/iT/i (1 ±^ex^) ■ (3.4) 

Clearly, this expression vanishes at A; = fc^ for x = 0, i. e., in the xy plane 
in /c-space. On the other hand, assuming 6e > 0, the |-Fermi sphere be- 
comes elongated along the z axis (^fc| vanishes at k > k^), while the |-sphere 



42 



Pairing with broken space symmetries: loff vs. dfs 



is squeezed. The conservation of the densities and pi requires a recal- 
culation of the chemical potentials, which is done by integrating again the 
corresponding momentum distributions. The net effect is that the surfaces 
approach each other on the xy plane, see Sect. I3.5]2l and Fig. 13.51 

3.3 The gap in the BCS, LOFF and DFS phases 

Consider a trap loaded with ^Li atoms and assume that the net number of 
atoms in the trap is fixed while the system is maintained at constant tem- 
perature. Assume further that the number of atoms corresponds to a Fermi 
energy Ep/kB = Tp = 942 nK, which in the uniforme and symmetric case 
at T = would translate into a density of the system p = 3.8 x 10^^ cm~^ 
and a Fermi momentum kp ~ 4.83 x 10^ cm~^. All the results below have 
been calculated for a homogeneous system at this density and at a constant 
temperature T = 10 nK <^ Tp, so the system is in the strongly degenerate 
regime. In the conditions of |MIT02| . this Fermi energy corresponds to about 
4 X 10^ atoms in a single hyperfine component of ^Li [see also |Had03| . espe- 
cially chapter 5]. Present experiments can control the partial densities in the 
two different hyperfine states |t) = |3/2,3/2) and |i) = |3/2, 1/2) of ^Li by 
transferring atoms from one to the other using ~76 MHz RF pulses |Had03] . 
Since the free-space triplet scattering length between ^Li atoms in these hy- 
perfine states is a = — 2160as, the system is in the weakly coupled regime 
kp\a\ <C 1. 

Without loss of generality, we assume a non-negative density asymmetry, 
i- e., Pi > Pi- We show only results for 6e > since we have checked 
that it is the one that gives the lowest free energy. We remind that this 
Se > corresponds to a cigar-like deformation of the t Fermi surface and a 
pancake-like deformation of the | Fermi surface, as we will see explicitely in 
Fig. 13.51 The pairing gaps of the LOFF and dfs phases computed from the 
(coupled) gap and number equations [ ()1.23|1 and (jSHJ] are shown in Fig. 13.11 

Table 3.1: The various candidates for a superfluid ground state studied. 



Candidate phase 


a P Se 


Symmetric BCS 
Asymmetric BCS 
LOFF 
DFS 



^00 
^0^0 
^0 ^0 
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Figure 3.1: Dependence of the pairing gaps in the loff phase (upper 
panel) and the DFS phase (lower panel) on the density-asymmetry parameter 
for the values of the the pair momentum P/kp and deformation parameter 
Se indicated in the legends. The Fermi momentum is kp = 4.83 x 10^ cm~^ 
and the scattering length is a = — 2160aB, so that kp\a\ = 0.55. 




as a function of the density-asymmetry parameter a = (p| — Pi)/{p] + pj,) 
for different values of the corresponding parameter signaling the symmetry 
breaking: total momentum P for the LOFF phase (top panel) and deformation 
5e for the DFS phase (lower panel). As expected, for vanishing asymmetries 
a 0, the maximal gap is attained by the standard BCS ground state, which 
is indicated by the solid, black line in both panels. The symmetry-breaking 
states that are very different from it (e. g., a loff state with P/kp > 0.1 or 
a DFS state with a > 0.1) have notably smaller gaps. 

The situation changes for increasing asymmetry. For example, the LOFF 
state with P/kp = 0.05 presents a larger gap than the asymmetric BCS 
state for a > 0.03, a density asymmetry for which also the DFS state with 
Se = 0.12 has a gap larger than the BCS one. Finally, for a > OmSf ~ 0.055, 
the BCS state no longer presents a finite gap, while both the LOFF and DFS 
phases 'survive' up to higher asymmetries: a^^^ ~ 0.06 (for P/kp = 0.05 — 
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0.10) and a^^^ ~ 0.075 (for 6e ~ 0.15), respectively. It is also remarkable 
the presence of a 'reentrance' phenomenon in the dps phase. For instance, 
consider the case 6e = 0.17 (blue, dotted curve). Such a large deformation 
shows a finite gap only for asymmetric systems, and a lower acri ~ 0.02 and 
higher acr2 ~ 0.07 critical asymmetries can be identified. As the quantity 
that determines the true ground state of the system is the free energy and 
not the size of the gap, we cannot tell from Fig. 13. II what will be the structure 
of the ground state. However, before discussing this point in detail, we shall 
study some aspects of the excitation spectrum of the system. 



3.4 Excitation spectra in the superfluid phases 



To elucidate the dominance of the phases with broken space symmetries over 
the asymmetric BCS state, it is useful to consider the modifications implied 
by these phases to the quasiparticle spectra 



These energies correspond to the poles of the propagators G-^i of chapter [H 
Physically, is the excitation energy of the system when we move it from 
its ground state to a state with momentum k and (pseudo)spin |. For a 
non-interacting system for which the number of particles is not considered 
fixed, this can be accomplished in two ways: 

• adding an j-particle with momentum fc, or 

• removing a |-particle with momentum —k. 

Conversely, —E^ is the excitation energy when the system is forced to have 
k momentum and | (pseudo)spin. Which is the process that ultimately will 
be required to actually perform these excitations is essentially determined 
by Pauli's principle and the interactions in the medium. For example, a 
system composed only of j-particles (/i = 5/i = /i|/2) at T = 0, can only 
be excited to states \k, t) for momenta > as for /c < all states are 
already occupied and Pauli's principle forbids a new j-particle to be added 
to the system below its Fermi momentum. The resulting system will have 
excitation energy h^k"^ /{2m) — ep = E^. For k < k^, the reachable states 
are \k, |) with excitation energy ep — h?k^/{2m) = —E^ > 0, corresponding 
to the creation of a hole in the |-Fermi sea at k. Similar ideas hold for the 
interacting system, and for other values of the densities p| and pj. 

We call the reader's attention to the minus sign in front of E^: this is due 
to the conventional way of assigning energies to hole excitations (see, e. g., 




(3.5) 
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|FW71| . pp. 74-75 ): the more negative they are, the more excited state of the 
system we get. This is easily understood again in the pure f-system, where 
removing an j-particle from k = leaves the system in its ground state, and 
therefore the excitation energy is eexc(^ = k^) = = —E^^. On the other 
hand, removing a particle at = will leave the system in a highly excited 
state, corresponding to having promoted [in the (A^ — l)-particle system] 
the t-particle at A; = to A; = the corresponding excitation energy is 
^excik = 0) = = . Therefore, the excitation energies of the system 
displayed in Fig. 13.21 are 

• eJ, := = -y^|~+~[Ap+ Ea- excitation energy of the system with 
momentum k and (pseudo)spin ]; 

• E^ := —E^ = Eg + I Ap — Ea- excitation energy of the system with 
momentum k and (pseudo)spin |. 

Let us turn now to the effects of the density asymmetry on the solution 
of the gap equation. We have seen in the previous chapter that, in the 
asymmetric BCS state, Ea acts in the gap equation ()1.23jl to reduce the 
phase-space coherence between the quasiparticles that pair. In other words, 
Ea 7^ introduces a 'forbidden region' for the momentum integration in the 
gap equation [cf. Sect. EIll especially Eq. ()2.2bll and Fig. [2121 . The BCS 
limit is recovered for Ea = 0, with equal occupations for both particles and 
perfectly matching | and | Fermi surfaces. This blocking effect is responsible 
for the reduction of the gap with increasing asymmetry and its disappearance 
above a = a^Sf — 0.055. 

When the pairs move with a finite total momentum or the Fermi surfaces 
are deformed (and taking the symmetry-breaking axis as z axis), the anti- 
symmetric part of the spectrum Ea is modulated with the cosine of the polar 
angle [cf. Eqs. ()3.3H3.4|) ]. In the plane-wave loff phase Ea oc x = cos 6*, 
while in the DFS phase, which is the object of our primary interest, we have 



This angular variation acts to restore the phase-space coherence for some 
values of x at the cost of even lesser (than in the BCS phase) coherence for the 
other directions. That is, the width of the forbidden region in Fig. 12.31 now 
depends on the direction: it is reduced in the xy plane and increased on the 
z axis. This effect can be explicitely seen in Figure ESI which compares the 
quasiparticle excitation spectra in the BCS and DFS phases. Let us comment 
carefully this figure, as it contain much information. The first column shows 



^5 = efc - At - 5iJ,5ex^ 
Ea = -Sn + udex^ ■ 



(3.6a) 
(3.6b) 
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the results for the usual symmetric BCS case; the second column has the 
results for the asymmetric BCS case with a moderate density asymmetry 
a = 0.04; finally, the third column contains the results for the DFS phase 
with the same asymmetry and the optimal deformation 6e = 0.08.* In all 
three columns, the top plot displays the excitation spectra i?^ (black lines) 
and eI (green lines) close to the Fermi momentum kp. Solid lines correspond 
to the results along the symmetry-breaking z axis, while the dashed lines in 
the last column stand for the spectra in the xy plane in /c-space. The figure in 
the top-left corner is equivalent to the lower panel in Fig. II. 21 The spectra for 
the asymmetric BCS case are shifted with respect to each other due to the fact 
that 7^ 0, cf. Eg. I3.5[ but keep the rotational invariance. Finally, the plot 
for the dfs phase shows the effect of the x-dependence of the spectra: the 
solid curves along the symmetry-breaking axis {i. e., x = ^ ov k± = 0) have 
gone further apart, while the dashed lines corresponding to x = ^ {kz = 0) 
have approached one another. 

To better evaluate the rotational properties of each phase, the lower plots 
show eJ} as a function of k^/kp and k^/kp. Keeping k = ^/k'^ + k'j_ constant 
amounts to moving along circles on the kz — k± plane, therefore exploring 
the angular dependence of the functions. For the symmetric BCS case, both 
I and t excitation spectra are equal, have a minimum at k = kp and are 
rotationally symmetric. For an asymmetric system (second column) the | and 
I spectra are no longer equally deep, but E^, is deeper than eJ,. Therefore, the 
quasiparticles of one species, that are defined near the corresponding Fermi 
surface, are far (in phase space) from those of the other species. The third 
column shows the excitation spectra in the DFS phase, as given by Eqs. 
and ()3.5|1 . Rotational symmetry is now broken, as the spectra along k^ (solid 
lines) are different from the spectra perpendicular to k^ (dashed lines). 

A most remarkable feature in the DFS is that the energy separation be- 
tween the quasiparticle spectra along the symmetry-breaking axis is consid- 
erably larger than in the asymmetric BCS state; in the orthogonal directions 
the opposite holds. Compared to the asymmetric BCS state, the phase-space 
overlap between pairs is accordingly decreased in the first region and in- 
creased in the second. The net result, displayed in Figure IT3l is the increase 
in the value of the critical asymmetry Omax at which superfiuidity vanishes. 
As noted above, at large asymmetries the DFS phase exhibits the re-entrance 
effect: pairing exists only for the deformed state between the lower and up- 
per critical deformations (ctcri 7^ 0). We note that to obtain this effect the 
recalculation of the chemical potentials through the normalization condition 
on the densities is essential, as it affects dramatically the value of 6fj,, that 



*By 'optimal' we mean 'with lowest free energy', see below. 
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Figure 3.2: Dependence of the quasiparticle spectra eJ, (black lines) and Ej, 
(green) on the momentum for the symmetric BCS case (left); an asymmetric 
BCS case with a = 0.04 (center); and the optimal 6e = 0.08 DFS result for 
the same asymmetry (right). The z axis is the symmetry-breaking axis of 
the DFS phase. In the top plots, the solid lines correspond to the spectra 
eJ, (black) and E^, (green) along the axis. The dashed lines represent the 
results on the xy plane {k^ = 0) with the same color coding. The dotted 
line is at i?fc = 0. In the lower plots, Ej, is given only for k^ > and E^, for 
< 0. 



BCS 



Asymmetric BCS 



DFS 




enters Ea- 

Notice also that in the asymmetric systems, the | spectrum is gapless in 
a region of momentum space defined by 

Ei = E^{k,,k±) <0 . (3.7) 

The possibility of exciting the system without energy cost has important con- 
sequences for the dynamical properties of the paired states, such as the trans- 
port and the collective modes, and leads to a number of peculiarities in their 
thermodynamics |Sa,r63l ISALOTf ISTTM ITWcH ISH03l ITTSMl IWY03| . That is, 
the macro-physical manifestations of the LOFF and DFS phases such as the 
response to density perturbations or electromagnetic probes and the thermo- 
dynamic functions (heat capacity, etc) will differ from the ordinary BCS phase 



48 



Pairing with broken space symmetries: loff vs. dfs 



Figure 3.3: Dependence of the critical asymmetry of the transition from 
the superfluid to the normal state on the pair momentum in the LOFF phase 
(solid line) and the deformation parameter in the DFS phase (dashed and 
dot-dashed lines). 



0.08 




0.05 0.10 0.15 

Pcr/kp. 6s 



due to the nodes and anisotropy of their spectra. We will show in Sect. 13.61 
how such an anisotropy can be used to discriminate phases with broken space 
symmetries in time-of-flight experiments. We remark finally that the phases 
with broken space symmetries (loff and DFS) present a larger number of 
excitation modes because of the breaking of global space symmetries; these 
additional modes are usually called Goldstone modes |Gol61[ lPS95| . 

3.5 Determining the ground state 
3.5.1 Calculation of the free energy 

The phase that must be identified as the equilibrium state at a given density 
asymmetry and temperature is the one that has the lowest free energy. We 
have calculated the free energy for each candidate phase (bcs, loff, dfs) 
in the following way. 

We have defined the free energy as usual by 

^ = -E'kin + Epot — TS, (3.8) 

where the first two terms comprise the internal energy which is the statistical 
average of the Hamiltonian ()1.4|) , T is the temperature and S the entropy 
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given for a gas of fermions by the well-known expression |FW71| . 

For a contact potential of strenght g , the potential energy is easily evaluated, 
and we have 

^kin + ^pot = j Cfc ht(fc) + n^{k)] - — , (3.9) 

n,(fc) = u\k)f{El) + t;2(fc) (1 - /(E-)) . (3.10) 

The free energy of the undeformed normal state follows by setting in the 
above expressions A = = 6e, while that of the BCS phase is given by 
A = Se, and so on for the loff and DFS phases (see Table EH}. 

Because of the contact form that we use for the interaction, the gap 
equation and the superfluid kinetic energy need a regularization. There are 
several ways to regularize the gap equation. We can write them formally 
together in the form [cf. Eq. (|1.23jl ] 

where fp is the usual Fermi distribution function. The case 7 = 1 and 
A — > cxD corresponds to the common practice of regularization |SHSH96] . 
which combines the gap equation with the T-matrix equation in free space, 
and is the one we have used in chapter 2. Choosing 7 = and a finite A 
corresponds to the cut-off regularization of the original gap equation. The 
appropriate cutt-off in this second scheme is found by requiring both schemes 
to give the same value for the gap. Then, this A is used to evaluate the kinetic 
energy contribution to the free energy. We note once more that Eq. ()3.1H1 
must be solved together with the normalization constraints on the densities 

/d'^k 
j^n.ik) (s=T,i). (3.12) 

3.5.2 Analytical estimation of the optimal deformation 

Before comenting the numerical results, it is instructive to perform an analyt- 
ical estimation of the deformation 6e of the Fermi surfaces that one expects 
will maximize the pairing energy. 

The physical idea behind the deformation of the Fermi surfaces of the 
pairing species is to approach them in some regions of momentum space, so 
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that the quasiparticles' phase space has a sizeable overlap. If this deformation 
had no kinetic energy cost, we can imagine that the optimal deformation 
would be the one that effectively makes the two Fermi surfaces be as close to 
each other as possible, given the constraint of the conservation of the different 
densities. Thus, we will calculate which deformation Se brings the outermost 
part of the Fermi surface of the least populated species and the innermost 
part of the Fermi surface of the majority to touch. The optimal deformation 
is somewhat different that the one predicted by this estimation because of 
the investment in kinetic energy that this deformation requires. In any case, 
we can get a good estimate of the optimal 6e in this simple way. 

We assume for simplicity T = and start from Eq. (j3.4|] . Here, we must 
understand fig = h^k1/{2m)^ where ks is such that the density of s-particles 
is conserved. For the case of j-particles: 

d^k 



1 ^1 

(27r)2 24 



/ —, ^ , arcsh v5e 

2Vl + 5e(5 + 26e) + 6- 



(3.13) 



The expression for Pi(fcj) is analogous with 6e —6e. As the 'true' Fermi 
momentum of j-particles is defined by p| = p(l + a)/2 = k'l/{6n'^), we have''' 



1 + a 

n 



k^ 1 — a 



_ k\ 2v/rT5i(5 + 25e) + 6^^^^^ ^k\ f 9 
~ kl 2^/^^(5 - 26e) + QSSS^M ^ kl\ 8 ^ 

where the last result follows for 5e <^ 1. In the plane perpendicular to the 
symmetry-breaking axis (cos^ = 0), each momentum distribution vanishes 
at the corresponding kg- Therefore, we will^have an exact matching of the 
two Fermi surfaces in this plane for k^ = ki- Then, for a given a we can 
obtain our estimation for the optimal deformation from the last equation: 



^Here it is clear that we are working at fixed density asymmetry, and not chemical 
potential difference, as was the original treatment of Larkin and Ovchinnikov |L0fi4| and 
Fulde and Ferrell JTfi^ for the loff phase, see below. In our approach, the chemical 
potentials come through the normaHzation of the densities, Eq. II3.12|I . which depends on 
the shape and location of the Fermi surfaces. 
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Substitution of this value into gives the new value for k^, and therefore 

the new chemical potential and analogously for fi^. For a -C 1 we expect 
a deformation Seopt ~ (16/9) a ~ 2a. We will see that this approximation 
agrees reasonably well with the numerical calculations below, even though 
the overlap is not perfect in the numerical solution because in this analysis 
we have disregarded the investment in kinetic energy. 

3.5.3 Numerical results 

We plot in Fig. 13.41 the difference AjF = JF5 — jF/y between the free energy 
J-'s of the superfluid phase (either BCS, loff or DFs) and the free energy jF/v 
of the normal state (A, 6e, P = 0). The BCS result is plotted in both panels 
as the black solid line, while the results for the LOFF phase are summarized 
in the top panel for various values of P/kp. The lower panel contains the 
corresponding results for the DFS phase for a variety of Je's. First of all, a 
simple comparison of this figure with Kill shows that A closely follows the 
behavior of — A^. However, the contribution from the kinetic energy to JF 
introduces some important differences as, for example, the delay in the DFS 
being preferable to the normal phase, e. g. for 5e = 0.17 from a ~ 0.018 
(where a non- vanishing gap appears) to a ~ 0.025 (where AjF < 0). 

The LOFF phase is preferred to the normal and the BCS phases in a narrow 
window of asymmetries 0.04 < a < 0.057 and for a total momentum of the 
pairs P/kp ~ 0.05, as shown by the red dashed line lying below the solid one. 
This is consistent with the results obtained in the scheme where the density 
asymmetry is described by fixing the difference in chemical potentials of the 
pairing species 6fi. In this approach, the critical value for the BCS phase is 
S^fcs ^ o.707A(0), while for the LOFF phase = 0.755A(0), where 

A(0) = A(5/i = 0) |TX)64I IFF64| . Note that while there is a non-trivial 
solution to the gap equation for P/kp > 0.1 the gain in pairing energy is 
less than the investment required in kinetic energy due to the motion of the 
condensate, and the net free energy of the LOFF phase is greater than that of 
the asymmetric BCS phase for these momenta. However the pairing energy of 
the LOFF phase could be improved by choosing a more complex form of the 
order parameter, e. g. by keeping a larger number of terms in its expansion 
in the Fourier series. 

The DFS phase is the ground state of the system in a wider range of 
asymmetries 0.03 < a < 0.06 for the deformation parameter in the range 
0.12 < 6e < 0.14. For even larger deformations the gain in pairing energy 
does not compensate the investment in kinetic energy due to the ellipsoidal 
stretching of the Fermi surfaces: the curve for Se = 0.17 lies always above 
the BCS one. We note finally that the free energy is also affected by the 
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Figure 3.4: Difference of the free-energy density of the plane-wave loff 
(upper panel) and dfs (lower panel) phases with respect to the normal phase, 
as a function of the asymmetry parameter a, for the values of the pair mo- 
mentum P/kp and the deformation parameter 6e indicated in the legends. 
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re-entrance effect («. e., restoration of pairing correlations as the asymme- 
try is increased), but in fact we do not find that this effect is responsible 
for the ground state at these weak couplings: the free energy for the case 
6e = 0.16 lies in our calculations always above those for Se < 0.15, which do 
not present re-entrance, in contrast to what happens in more strongly-coupled 
systems such as neutron-proton pairing at saturation density |MS 02] or 'two 
color superconductivity' between up and down quarks in dense hadronic mat- 
ter !MS03a| . 

To summarize, the coherence is restored and the strength of pair-corre- 
lations is increased in the loff phase due to the finite momentum of the 
Cooper pairs. In the dfs phase, the same is achieved by stretching the 
spherical Fermi surfaces into ellipsoids. The fundamental difference between 
these phases is that translational symmetry remains intact for the DFS phase, 
as it breaks only the rotational symmetry, while the LOFF phase breaks both 
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symmetries. Quantitatively, the maximal value of the gap and the absolute 
value of the ground-state free energy are larger in the DFS phase than in the 
LOFF phase for asymmetries a > 0.04. For these asymmetries both phases 
are favorable over the homogeneous BCS phase. However, one should keep in 
mind that the LOFF phase admits a variety of lattice forms, and the plane- 
wave structure might not be the most favored one |ABR,OH ISedOU IBR,02[ 

3.6 Detecting the DFS phase in experiments 

Experimental evidence for the phases with broken space symmetries can 
be obtained from the studies of their momentum distributions which, un- 
like in the homogeneous phase, must be anisotropic in space. Indeed, from 
Eqs. ()3.5p . ()3.6p and ()3.10|l . we see that, for 5e ^ 0, the probability of a 
s-particle having momentum fc, ns{k), will depend on k"^ and k\ separately 
or, equivalently, on k = \k\ and x = cos9 = kz/k, and not only on the scalar 
combination k"^ = k"^ + k\. 

In this case, it is useful to analyze the angular anisotropy in momentum 
space of the occupation numbers, that we define as 

5n,{k) = Mk; X = 1) - n,{k; X = 0)| . (3.15) 

This quantity vanishes for rotationally-invariant systems, and can be taken as 
a measure of the degree of symmetry breaking of a given phase. For example, 
5ns{k) = (Wk) for the symmetric and asymmetric BCS solutions to the gap 
equation. For the DFS phase, the momentum distributions for each species 
depend on x, and we expect that 6ns will not vanish but present a maximum 
around the corresponding Fermi momentum kg. 

We show in Fig. 13.51 the momentum distributions for the same system 
as in Fig. 13.11 kp = 4.83 x 10^ cm~^ (indicated by the vertical line) and 
a = — 2160as. For clarity, we present only the data around kp, where all the 
interesting physics happens. There are four sets of curves corresponding to 
the symmetric BCS phase (thin solid line), asymmetric BCS with a = 0.04 
(dashed lines) and DFS with the same asymmetry and 6e = 0.08 for x = 
(dash-dotted lines) and x = 1 (dotted lines), respectively. For the a ^ 
cases, blue lines correspond to ^i(fc) while red lines are for ni{k). Finally, 
the thick, bell-shaped curves correspond to Susik) with the same color coding. 

The Fermi surfaces of the symmetric BCS phase are spread around kp over 
a width of the order of kp A/ep. One should notice that some broadening of 
the surface is also due to the finite temperature assumed in the calculations. 
The asymmetric BCS distributions (dashed lines) show a similar shape, but 
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with the characteristic decays centered around the corresponding Fermi mo- 
menta ki ~ 4.75 X 10^ cm~^, fc^ ^ 4.91 x 10^ cm~^. At the low temperature 
studied, the momentum distributions still fall off very fast, and the two Fermi 
surfaces are separated ~ 0.16 x 10^ cm~^, thus triggering the Pauli blocking 
effect commented in Sect. 12.21 

The situation changes when looking at the DFS results which show a clear 
dependence on x- First of all, we notice that the DFS curves for x = (dash- 
dotted lines) are almost indistinguishable, while the results for x = 1 (dotted 
lines) are further apart. This means that the (ellipsoidal) Fermi surfaces of 
t and I particles have approached each other in the xy plane in momentum 
space, while they have gone apart along the z direction (see Fig. 13 .511 . Note 
also that the range of occupations covered by the DFS curves includes the 
asymmetric BCS results. One could have expected the asymmetric BCS results 
to coincide with the DFS ones for x = 0, as in both cases the symmetry- 
breaking term 5ex^ vanishes [see Eq. ()3.6I1 ]. However, the fully self-consistent 
solution of the gap equation together with the number equations involves a 
re-calculation of the cheniical potentials for the DFS case, resulting in new 
effective Fermi momenta kg (see Sect. 13. 5. 2^ which become very similar to 
each other and to kp- That is the reason why the DFS results at x = are 
so close to each other, and relatively far from the asymmetric BCS ones. In 
addition, we obtain a larger gap (cf. Fig I3.1jl . that translates into slightly 
softer decays of the DFS momentum distributions as compared with those for 
the asymmetric BCS case. 

The anisotropy in the DFS occupation probabilities along the directions 
parallel and orthogonal to the symmetry-breaking axis reaches a maximum 
of more than 90% for both species centered at different momenta: for Sn^, 
around k^y/l — 6e ~ 4.63 x 10"^ cm~^ [momentum for which 5nj(x = 0) ~ 1 
and Sni{x = 1) ^ 0]; for at k^^/T+~6e ^ 5.07 x 10"^ cm~^ [where 
^^t(x = 0) and 5n|(x = 1) «i 1]. 

In practice, this means that one expects to find more t particles with 
momenta k > kp in the z direction (x = 1) than in the radial direction (x = 
0), while the opposite holds for the | particles. Thus, a direct way to detect 
the DFS phase is a measurement of this crossed anisotropy in the momentum 
distributions of the trapped atoms. Such a measurement can be realized by 
the time-of-fiight technique |.TTT; A95L ImTT951 IRice99] . This method uses the 
fact that after switching off the trap, the atoms fly out freely and an image of 
their spatial distribution taken after some time of flight provides information 
on their momentum distribution when conflned inside the trap. Assuming 
that the system was in the deformed superfluid state one would detect a 
number of particles of type t (majority) in the direction of symmetry breaking 
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Figure 3.5: Dependence of the occupation probabilities of two hyperfine 
states on the momentum. The Fermi momentum kp = 4.83 x 10'^ cm~^ 
is indicated by the vertical line. The labeling of the lines is as follows: 
a = = 6e (black, solid line), a = 0.04 and 6e = (dashed lines); a = 0.04, 
6e = 0.08, X = (dash-dotted) and x = 1 (dotted). The blue lines are 
for I particles (major population, see left), while the red ones are for | 
particles (minority). The bell-shaped curves show the anisotropy as defined 
by Eq. (l3T5ll for a = 0.04, Se = 0.08. 
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by about 90% larger than in any orthogonal direction, and viceversa for the 
I particles (minority). Therefore, the presence of this crossed anisotropy in 
the detected momentum distributions would be an evidence for a 'deformed 
superfluid state' being the ground state of the system, as deformation alone («. 
e., without pairing) would not lower the energy so as to produce a deformed 
non-superfluid ground state. Note that this argument is equally valid for 
a homogeneous system or for an atomic sample in a spherical trap, where 
no privileged direction is introduced by the trapping potential. For a non- 
spherical trap, the normal-state momentum distributions of both species are 
expected to be deformed in the same way. Therefore, the detection of a 
crossed anisotropy in the momentum distributions as discussed above would 
also be a strong case for the DFS phase being the ground state of the system. 
However, a more specific calculation (e. g.., in local density approximation in 
the trap) would be necessary to quantify this effect and study the influence of 
the trap anisotropy on the momentum distributions in the different phases. 

We remark finally that the direction of spontaneous symmetry break- 
ing (in fc-space and, therefore, also in real space) is chosen by the system 
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randomly and needs to be located in an experiment to obtain maximum 
anisotropy. Also, a clear distinction between the dfs and the loff phases 
can be achieved in the time-of-flight experiments, since the latter predicts 
periodic momentum distributions. 

3.7 Summary 

In this chapter we have studied the possibility to generalize the BCS ground 
state by letting the Cooper pairs carry a finite center-of-mass momentum 
(the so-called loff phase) or by deforming the Fermi surfaces of the pairing 
species (dfs phase). For the density-symmetric case (a = 0), the BCS pair 
wave-function is the best among these three options, i. e., it gives the lowest 
free energy for the system. However, for systems with different populations 
of the two species (a 0), this is only true for very small asymmetries 
(a < 0.03 = 3%). For larger asymmetries, the loff and/or dfs phases be- 
come preferable. Our quantitative analysis has shown that, at weak-coupling, 
the DFS phase is prefered to the normal, BCS and plane-wave loff phases 
(though more general forms for the LOFF phase are possible and might lower 
further the free energy). We note that similar conclusions have been ear- 
lier reported for strongly-coupled systems such as neutron-proton pairing in 
asymmetric nuclear matter in the ^Si—^Di channel |MS02L[MS03b| and pair- 
ing of up and down quarks of two colors in dense hadronic matter |MS0 3a,]. 
Therefore, we expect a similar behavior for more strongly-coupled atomic 
systems, such as those produced by means of magnetic Feshbach resonances 
to stu dy the BCS-BEC crossover |.TTT. A03a[ IMFTM IRiceOSi ITnns04l IENS04al 
iDui^OBj. 

We have finally shown how the dfs phase can be detected in an ultra- 
cold gas of fermionic atoms in a spherical trap by studying the momentum 
distribution of the released atoms in a time-of-fiight experiment. 



Chapter 4 

Pairing in boson-fermion mixtures 



His house was perfect whether you liked food, or sleep, or work, 
or story-telling, or singing, or just sitting and thinking best, or a 
pleasant nnixture of thenn all. 

J. R. R. Tolkien, The Hobbit 

4.1 Introduction 

Mixtures of quantum fluids of different statistics have been an interesting 
field of research for a long time, specially in the context of helium systems, 
where the fraction of ^He (fermion) in a homogeneous ^He (boson) medium 
is limited to a;3max ~ 6.5%. For X3 > a;3max the system becomes unstable and 
separates into a mixed phase with a ^He concentration xamax and another 
phase with pure ^He atoms. 

Since the first achievement of Bose-Einstein condensation, the purpose of 
building an analogous ultracold Fermi system had to overcome the cooling 
limitations established by Pauli's exclusion principle, that forbids s-wave col- 
lisions between indistinguishable fermions. This problem can be solved in a 
boson-fermion mixture, as the bosonic component can be 'easily' cooled and, 
by thermal contact, drives the fermionic component down to ultralow temper- 
atures, as was first achieved at Rice University when a gas of ^Li was driven 
to the degenerate regime by sympathetically cooling it with ''Li jRiceOlj . 

After this experiment, many others have shown the feasibility of study- 
ing interesting quantum phenomena in ultracold fermionic systems. For ex- 
ample, the LENS group used ^^Rb-^°K mixtures with varying numbers of 
atoms of both species to study the collapse predicted by mean-field the- 



58 



Pairing in boson-fermion mixtures 



ory IM01981 lR.Fn2] and extracted a value for the interspecies scattering length 
a(^^Rb-^°K) = -22 nm [LENS02J.* 

In this chapter, we address the problem of determining the fermion- 
fermion interaction inside a dilute mixture of bosons and fermions. This 
problem is important as medium effects can have dramatic consequences on 
the behavior of the system. A classical example is that of superconductivity 
in metals, where the phonon-mediated interactions between electrons give 
rise to an attraction between the latter, thus triggering the Cooper instabil- 
ity of the system and resulting in the superfluid behavior of the electronic 
component [Sch88j. First, we will review the three-dimensional situation, 
already studied by Viverit et fl/. |VPSnni IRHSnnj . Then, we will study the 
two-dimensional case, where the reduction in the dimensionality of the phase 
space makes one expect that correlations will play a more important role 
as compared to the 3D situation. Indeed, today's available capabilities to 
manipulate trapped, atomic gases allow for the production of effectively one- 
and two-dimensional systems by the deformation of the traps or, more effi- 
ciently, by shining on it an array of standing light- waves that produce a kind 
of periodic potential known as optical lattice (see, |GV99| INISTn2| and ref- 
erences therein). We will end up applying the developed model to the cases 
of ^^Rb-^°K and ^Li-^Li mixtures under current experimental investigation. 



4.2 Three-dimensional mixtures 

Let us start by considering the fermion-fermion {FF) interaction in three- 
dimensional space for a one-component fermionic system in the presence of 
bosons. At the ultralow temperatures of interest T < 10/iK and dilute con- 
ditions of present experiments, only s-wave collisions are relevant (unless 
they are absent). As for a system of indistinguishable fermions s-wave in- 
teractions are forbidden by Pauli's principle, we will be left only with the 
boson-fermion {BF) and boson-boson [BB] interactions, while the direct 
FF one can be neglected. For a dilute and ultracold system, these collisions 
can be modeled in three dimensions by a contact interaction or, equivalently. 



*This value is currently under revision as determinations by similar, as well as other 
presumably more precise, methods have given a more moderate value a — —15.0 ± 0.8 
nm j ,TILAn4allBon05] . 
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by momentum-independent T-matrices, 

Tbf — dBF 

rriBF 

_ inh^ 

Tbb — ciBB , 

niB 

where (mp) is the mass of a boson (fermion) and uibf = niBmp / {mB + 
mp) is the reduced mass in a boson-fermion collision. 

The boson-induced fermion-fermion interaction can be represented by the 
diagram in Fig. l2.2r bV where now the bubble corresponds to a density fluc- 
tuation in the bosonic medium. Analogously, we can use a formula similar 
to Eq. (lOl) : 

(fcTF|fc) = n|p^(|fc'-fc|)T|p. (4.1) 

Here 11^^^ stands for the density-density response function of the bosonic 
component in the random phase approximation (RPA). Already from its 
diagrammatic representation [Fig. I4.2f b)]. one can see that it is possible to 
evaluate n^^"^ as a series once Tbb and the bosonic Lindhard function lis are 
known. For in- and out-going fermions on their Fermi surface there will be 
no net energy transfer through the bosonic medium [Fig. l4.2r a)]. and one can 
take the static limit a; = for the bosonic response function [NPQO, .FW71| . 

nB(g;u; = 0) = ^-^ . 



Altogether, we get 



Note that, due to the low density of the systems we are interested in, we are 
neglecting any fermionic influence on the bosonic component. 

As in Sect. 12. 4L we must now project this polarization-induced interaction 
onto its p-wave contribution with the corresponding Legendre polynomial. 
Again, for a weakly-interacting system it is sufficient to consider only the 
case where both in- and out-going fermions are on their Fermi surface: 

V^^='\kp, kp) = T|p / -zUY'^{yW^)kF) = -^hj,{x) , (4.3) 

J-l ^ I BB 

with X = li^kp/ {niBPB'TBB) = '^{^ f / PB){jnF I f^B) [for a dilute, homogeneous 
bosonic system, hb = TbbPb], and the function 
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Figure 4.1: The functions governing the boson-induced fermion-fermion 
interaction in three (lower curve) and two (upper curve) dimensions. The 
dashed lines indicate the position and value of the corresponding maxima. 
The inset shows h2T) together with the quadratic approximation around its 
maximum (thin line with crosses). 
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which is plotted in Fig. 14.11 

This function presents a broad maximum at (x^^ ~ 3.48, /i^^ ~ 0.10). 
This means that the induced interaction is optimized for a particular bo- 
son/fermion ratio given by 

2/3 

^ 2.88 asB , 

PB 

or, in terms of energies, 

^^1.74^. 

For example, for the case of mixtures of ^^Rb (boson) and '^°K (fermion) 
(o-BB = 5.2 nm |vKKHV02] ). this relationship gives pp = 1.8 x 10^^ cm~'^ 
for the case ps = 10^^ cm~^, or ep = 138 nK for pB = 36 nK. Estimating 
the resulting p-wave gap for the optimal boson concentration by means of 
Eq. ()2.9|1 . we get for these densities Ai/ pp ~ 10"'', which corresponds to an 
unattainable critical temperature. 
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4.3 Two-dimensional mixtures 

Pairing in two dimensions has the pecuhar feature that, for an attractive 
s-wave interaction between two different fermionic species, a two-particle 
bound state (of binding energy Eb) is always present, no matter how weak 
the attraction is. Therefore the system enters the strong-coupling regime 
even at low densities [SRVR89, RDSHHI IH?l)S90t IPBS03| . forming a Bose 
condensate of fermion pairs characterized by 

Hf ^ £f — Eb/2 , 



A, ^ ^/2EbeF , 

where is the s-wave pairing gap in 2D. This is in clear contrast to the three- 
dimensional case for two reasons: (a) in 3D, a two-body (i. e., in a vacuum) 
bound state requires a minimum strength of the (attractive) potential, and 
(b) a bound state in the many-body system is only very weakly bound in the 
low density limit (the gap vanishes exponentially with kp 0). 

In three dimensions a small excess of one type of fermions implied an 
important reduction of the gap, and eventually its disappearance even for 
very small asymmetries (see Sect. I2.2ll . The situation in the strong-coupling 
regime is very different. Here, the system will form all possible pairs, while 
the remaining particles just stay in their original unpaired states, in close 
analogy with the behavior of nuclear matter at densities low enough to allow 
for the formation of deuterons (deeply bound neutron-proton pairs) |LNS"'"fll| . 

We therefore exclude in the following this 'trivial' case, and focus on 
treating a system of identical (spin-polarized) fermions. The first possibility 
of pairing concerns the p-wave pairing gap, Ai = AL^i^kp), which in the 
low-density limit is given by the weak-coupling result jRDSQflj 



Ai 

— = ci exp 

fJ,F 



27r^2 



where Ci is a constant of order unity and 



(4.4) 



TF = T^F~'\kF,kF-2fiF)= r^cos0(fc'|T(2/i^)|fc) , (4.5) 

Jo 

|A;| = |A;'| = kp , cos(j) = k' ■ k , 

is the relevant T-matrix element of the interaction, computed at collision 
energy 2/xj?. Note that the two-dimensionality of the system is reflected 
already in the integration to be performed over the angle 0, in contrast to 
the 3D case, where we could integrate directly over z = cos(0) [cf. Eq. ()4.3|) ]. 
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Figure 4.2: (a) Polarization interaction T between two fermions (dashed 
lines) mediated by the presence of bosons (solid lines). The labels indicate 
the momentum and energy of each line. For condensate bosons and fermions 
on the Fermi surface, ^ = 0,ci; = 0. (b) Diagrams contributing to the boson 
bubble in RPA; the last one is an example of a backwardgoing diagram, 
negligible when /i^ ^ 0. Here, thick solid lines are full propagators, thin 
solid lines are free propagators, and wiggles represent interactions. 
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Let us consider a homogeneous mixture of bosons and fermions, and dis- 
regard for the moment the possibility of a direct p-wave interaction between 
the latter. The most important contribution to the FF interaction is then 
given by the density fluctuations in the boson condensate sketched in Fig. 14. 21 



4.3.1 Constant T-matrices 

As a flrst approach to the problem, we take Tbf and Tbb to be constants as 
in three dimensions. Thus, the treatment is completely analogous to that of 
Sect. 14. 2| and we just have to evaluate the integral in Eq. (j4.5jl . We obtain 

^tt = T|f r - COS0 n*^(v/2(l-cos0)M = -^hj, (x) , (4.6) 

Jo I BB 

where x = h'^kp/ {itibTbbPb) as before, and 

^ / \ 2 I X 2 , , 

/i2D X := ^== - - . 4.7 

This function has a maximum located at (x^p^ = 2/(^/2 — 1) ^ 4.828, /i^pt = 
3 — 2\/2 0.172), and can in its vicinity be approximated by a parabola, as 
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shown in the inset of Fig. 14.11 This translates into a sharp Gaussian peak 
for the gap function, according to the weak-coupling result ()4.4jl . 

Note that, even though and are both dimensionless and simi- 
lar in magnitude, the change in dimensionality translates into different re- 
lationships of these parameters with the densities. In fact, x oc p'l^ j Pb in 
three dimensions while x oc Pf/ pB in 2D. Thus, the maximum boson-induced 
fermion-fermion attraction would be reached in two dimensions for a density 
ratio 

^ = ^ 2 ^ 0.7685^^ (4.8) 

if the collision T-matrices were considered constant. 



4.3.2 Energy-dependent T-matrices 

In two dimensions the s-wave scattering matrices Tbf and Tbb are not con- 
stant for a vanishing collision energy, but present a logarithmic dependence 
on the c.m.s. energy E of the two-particle state [ Arlh86i [M77Rn2] . i. e., 

(mPcM^O.B^0)|*)^^j-^. ,4.9) 

where rrired is the reduced mass of the colliding particles and i?o ^ is a 
parameter (with dimensions of energy) characterizing low-energy scattering. 
Therefore, we need to determine what are the possible collisions occurring 
in the mixture and what is the scattering energy corresponding to each one. 
As we consider elastic collisions — i. e. no change in the internal state of the 
atoms is permitted — , the relevant energy for the collision will be given by 
the square root of the invariant E^^^ = Q^Q^ {Q is the total 4-momentum of 
the collision), minus the rest energy. 

As indicated above, we are initially neglecting direct FF interactions, so 
we need to consider only the following events (shown diagrammatically in 
Fig.Hni): 

(a) a fermion on its Fermi surface (i.e., with 4-momentum [hk^Ep = 
\/vr?p + W'k'^p) in the laboratory frame) scatters off a boson in the con- 
densate (with 4-momentum (0,771^)). 

The total momentum is Q = {hh, Ep + m^), and so, assuming hkp <^ 
mp, we obtain 

1 ruB h'^k% ruBP 

-Einv = niB + nip + ■ =^ Esc = £f , (4.10a) 

ZTTiB + mp mp mp 
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Figure 4.3: Possible collision events in the mixture, according to Eq. (j4.1()jl . 
Dashed lines denote fermions, solid lines bosons, and wiggles represent inter- 
actions. The labels indicate the momentum and energy of each particle. 

k,eF k\eF 0,0 q,0 q,0 0,0 

0,0 q,0 0,0 -q,0 0,0 q,0 

(a) (b) (c) 



with rriBF = rnBrnp / {jriB + rnp) the reduced mass of the boson-fermion 
system. 

(b) two condensate bosons {Q^uib) scatter off each other. Here the calcu- 
lation is simple: 

Esc = . (4.10b) 

Because of this, the contribution of the backwardgoing diagrams like 
the last one in Fig. I4.2r b) vanishes. 

(c) a boson in the condensate scatters off a non-condensate boson with 
4-momentum {hq,EB = ^/rn^~+^?(f) . After removing a term 2mB, 
we get 

E,, = -^ {hq^ms). (4.10c) 
4ms 

Thus, the transition matrices now depend explicitly on energy: 

^ , , Anh^ 1 

Tbb g = 7-7-. , 4.11a 

rriB \n{AmBEBB/ri^q^) 

J {k )- ^^^^ ^ (4 lib) 

where Ebb and Ebf are the parameters characterizing low-energy s-wave 
BB and BE collisions, respectively. 

As noted above, within the approximation fiB = (or more precisely 
I^b -C fJ,F), events of type (b) do not contribute [Tbb(-£' = 0) = 0] and 
should therefore be removed from the series defining n^^"^. In practice, 
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this is performed by replacing TbbHb TBBnB/2 in the denominator of 
Eq. (I4.2jl . Thus, the polarization-induced interaction reads 

Tr'\kFM) = -"^^^^^i^hixX) , (4.12a) 
h{x,C)= 1 (4.12b) 

Jo TT X(1-COS0) - 

where we introduced x = Ji^kp/ (Atipb) = pf/pb, and the following dimen- 
sionless parameters 



CikF) 



h kpniBF mBF £f 



2'mp rriF mp Ebf 



f}^k% 1 mp Sp 



2mF Ebb ttlb Ebb 



These give a measure of the characteristic energy of the fermions relative 
to the scattering parameters Ebf and Ebb- We remark that the condition 
C ^ 1 must hold for the use of Eq. ()4.9p to be valid. 

Given a fermion density (and, therefore, C), the maximum value of h is 
reached for an optimal ratio a;opt(C)- For C 0, one obtains the following 
quasi-linear dependences on ln(^: 



X 



opt 



Ihj) -0.58693 [inC - 0.35197] , (4.13a) 



KAy) ^ 0.52022 [inC - 4.3257] . (4.13b) 

The accuracy of these approximations can be observed in Fig. l4.4l that shows 
the numerically computed optimal value Xopt(C) [Isft panel] and the corre- 
sponding /iopt(C) = ^(^^opt^C) [right panel] as a function of In^ (solid lines), 
compared with the approximations above (dashed lines). One readily notices 
the very good agreement between both calculations, especially for /lopt, for 
which already at C ~ exp(— 2) ?a 0.13 the two curves become indistinguish- 
able at this scale. 

Taking all these facts into account, the value of the pairing gap under 
optimal conditions becomes 

\}^ (4.14) 

CiPf mBmp 0.52022 [in C - 4.3257] ^ ' 



This means that the position and value of the maximum induced interaction 
will depend logarithmically on the Fermi momentum. 
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Figure 4.4: The optimal values Xopt and /iopt for the pairing interaction, 
Eq. (|4.12|) . The dashed lines indicate the asymptotic behaviors, Eq. ()4.13p . 




In C In C 



We compare now the importance of this induced interaction with a pos- 
sible direct fermion-fermion p-wave interaction, which at low density is given 
by |RI)S90| . 

(L=l) , Ah"^ Hf 

TJlF -C/1 

where Ei is the parameter characterizing 2D low-density p-wave scattering. 
As the boson mediated attraction, Eqs. ()4.12I4.13|) . at low density depends 
only logarithmically on the fermion density, it will dominate in this limit. For 
the same reason, any fermionic polarization corrections can also be neglected. 

Finally, we analyze the assumption ps <^ /Uf used above. The boson 
chemical potential is determined by |Sch71[ IFH88. .SR93t IT:MDR02| 

Anh^PB 1 

Pb = Pb \ bb{^ = Pb) ~- 



rriB In {Ebb /pb) 
while for free fermions we have 

27rh'^PF 



Pf 



Since the logarithm in the low-density domain is always large, we have the 
sufficient condition 

2mF pB ^ , Pf ^ ^uif 



TUB Pf Pb ^b 



which will be well satisfied in atomic mixtures {2mF/mB ~ 1) in the regime 
of validity of ipTl . 



4.3 Two-dimensional mixtures 
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Figure 4.5: The pairing gap for optimal boson concentration, Eq. (I4.14|) . 
as a function of i] = uibf^f I {jtifEbf) and C = t^f^f I iimBEBB)'- (left) for 
a ^^Rb-^°K mixture; (right) comparison with a ^Li-^Li mixture for fixed C,. 




4.3.3 Prospects of experimental detection 

In order to estimate typical sizes of the expected gap, we plot in Fig. 14.51 the 
gap Ai/^p^ according to Eq. ()4.14p . as a function of the parameters r] and 
( defined above (assuming for simplicity Ci = 1). On the left we show the 
calculation for a ^^Rb-'^^K mixture. We see that quite large gaps Ai ~ /i^? 
are achievable, and that it is mainly t] oc Sf/Ebf that determines the size 
of the gap. Indeed, the lines of constant value of the gap are almost parallel 
to the C axis, as shown by the contour lines in the r/ — C plane, except for 
Tj > 0.1, when the gap is already large {Ai/ fiF > 0.8). On the right part 
of the figure, we compare these results for the ^^Rb-'^^K mixture with the 
corresponding ones for a ''Li-^Li mixture. The general behavior is essentially 
the same, but the gap is a little bit smaller in the last case, the difference 
being due to the different value of the masses' coefficient in Eq. ()4.14p . whose 
value is 0.22 for the first mixture and 0.25 for the second one. Even though 
this difference is small, the fact that appears in the exponential gives rise to 
the difference in the figure. In summary, we expect quite large gaps Ai < fiF 
with fermion chemical potentials /iF ^ Ebf, while Ebb should play a minor 
role. 

In order to translate this condition into experimental quantities, we use 
the results of Refs. [PHSOOl ILM l)H02| . relating the 2D scattering pa- 
rameter Eq to the value of the 3D scattering length a^o, for a system confined 
in a strongly anisotropic trap characterized by frequencies ujj_ and uj^ ^ uj±. 
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Assuming that these frequencies are the same for the two trapped species, the 
treatment of |PS01| is equally valid for a mixture if we substitute the masses 
m there by twice the reduced mass of the pair. Thus, for the boson-fermion 
collision parameter, we obtain 



where B ^ 0.915 and Iz = h/ \/2mBF^z- Since at the same time for a 2D 
situation the condition jjLp <^ Huz must be fullfilled, one can only expect 
observable gaps if the exponential term is not too small. We can distinguish 
two cases: 

(i) aBF > 0: in this case the ratio Iz/clbf should be minimized as much 
as possible. One possibility is to strongly compress the trap in the z 
direction, or even better use a one-dimensional optical lattice to divide 
the trapped gas in a set of quasi-2D sub-systems. Another option 
is using a Feshbach resonance to enhance the repulsion; in this case, 
problems of stability might appear, as for strong interspecies repulsions 
the system may prefer to spatially separate the two species. 

(ii) asF < 0: in this case the exponential term is never small and pair- 
ing can be expected as long as fip/i^z) is not too small. Using the 
Thomas-Fermi approximation fiF = \/2NFhuJ±_ for the chemical po- 
tential of a non-interacting, strictly two-dimensional Fermi gas in a 
(in-plane) harmonic trap of frequency cj^, this last condition can be 
expressed as 



Let us take typical values from |LENSfl3] . where around 5x 10"^ atoms of 
"^^K were sympathetically cooled with ^^Rb down to T/Tp = 0.3 (Tp = 
430 nK) in a magnetic trap with frequencies u;^ = 27r x 317 Hz, Ua = 
271 X 24 Hz, and then submitted to an optical lattice with associated 
harmonic frequencies in the minima uji = 2t[ x 43 kHz. Therefore 
<^±/<^2 ^z/^i ~ 1/136, and we get a value ~ 2 for the ratio in 
Eq. ()4.16|) . which is close to the requirement for 2D behavior. 




(4.15) 



hhJz 



^/2N~F- 



(4.16) 



Let us note also that Eq. ()4.15jl applied to quasi-2D, bosonic ^^Rb trapped 
with the same Uz as above gives hfIEbb oc ^ ^ 1, thus ensuring a broad 
range of hfI Ebf over which large gaps are expected (see Fig. I4.5|l . 
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4.4 Discussion 

Under favorable circumstances, p-wave pairing gaps of the order of the Fermi 
energy seem to be achievable, and comparable to those predicted for s-wave 
pairing in quasi-2D two-component Fermi gases |PRSn3| . However, more 
precise quantitative predictions cannot be made in this regime within the 
present approach since, when iip ^ {Ebf, Ebb}, also the asymptotic ex- 
pression Eq. ()4.9|) becomes invalid. We stress that the same effect in three 
dimensions is less efficient in increasing the size of the gap, and one expects 

A,/fiF < 0.1 |HPSvnn| . 

Also the stability of a Bose-Fermi mixture in two dimensions is as yet 
an unexplored subject. In three dimensions, this topic has been studied by 
differe nt authors |M0l98, VPSnfH iMSYfVll \RFM rrZWMn2[ ISPZMlM IR otn2[ 
IVGn2j that reach similar conclusions. Here we briefly comment on the im- 
plications from these studies that might be applied to our case, but a more 
precise analysis of the two-dimensional case would be of high interest. 

According to Ref. |VPSOO| . in a three-dimensional boson-fermion mixture 
one can expect to find one of the following situations: 

(i) a fermionic phase and a bosonic phase, separated from each other; 

(ii) a fermionic phase and a boson-fermion mixture; 

(iii) a single uniform mixture. 

In case (i) there is no boson-fermion induced interaction nor sympathetic 
cooling. In case (ii) these problems are overcome, but only a fraction of 
the fermions is efficiently cooled and can undergo the superfluid transition. 
Therefore, the interesting situation is that of case (iii). This can be obtained 
if there is attraction between bosons and fermions (to avoid their spatial sep- 
aration), but in this case the system may collapse due to this same attraction 
as predicted in |RFn2j and observed in |LENSn2] . This will happen if, e. g., 
the number of bosons exceeds some critical number N^r, which will depend 
on ttBB and Gbf- For a uniform system, we know that Qbb > is required in 
order to avoid the collapse of the bosonic component. Roth and Feldmeier 
have shown that this condition also stabilizes significantly the mixtures, even 
for qbf < 0, while the case o^f > rapidly gives rise to spatial separation 
of the two components |RFn2j . 

Applying these arguments to the mixtures used in typical atomic ex- 
periments, we see that the case ''Li-^Li with gbb = —1-5 nm and gbf = 
2.2 nm |Riceni| does not correspond to the optimal stability conditions. 
However, the presence of the trapping stabilizes the system so that exper- 
iments can be performed. On the other hand, for the ^^Rb-^°K mixture. 
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where gbb = 5.2 nm |vKKHV02| and a^F = -22 nm |LENS02| . the stability 
conditions for the homogeneous case are fully satisfied. 

In conclusion, the ^^Rb-^°K mixture seems to be the best available can- 
didate within present-day ultracold boson-fermion mixtures to explore the 
outlined possibility of a p-wave superfluid transition in a (quasi)-2D system, 
both because of its demonstrated stability against phase separation and col- 
lapse, and because the expected gaps Ai ~ fip are larger that those for 
^Li-^Li mixtures. 



Chapter 5 



Dynamics of spin-1 condensates 
at finite temperatures 

Fu mio padre il primo ad accorgersi che qualcosa stava cambiando. 
lo ero appisolato e il suo grido mi sveglio: 
— Attenzione! Qui si tocca! 

Sotto di noi la materia della nebula, da fluida che era sempre stata, 
cominciava a condensarsi. 

Italo Calvino, Sul far del giorno (Le Cosmicomiche) 

5.1 Introduction 

In the previous chapters we have studied the prospects for a superfluid tran- 
sition in a low density fermionic gas with two species of fermions (chapters El 
and EI) or in a mixture of a bosonic and a fermionic species (chapter lij. We 
will analyze now the behavior of another multicomponent system, this one 
with only bosonic components, namely a so-called spinor condensate. 

We have seen in Section EIU how alkali atoms can be trapped by means of 
an inhomogeneous magnetic field. In this situation, the same magnetic field 
that traps the atoms by coupling to their magnetic moment, freezes their spin 
degree of freedom. Thus, the atoms behave in practice as spinless bosons, 
and the structure and dynamics of low density gases can be well described by 
the famous Gross-Pitaevskii equation for a scalar order parameter ip{r, t) = 
(\l/(r,t)) in an external potential Kxt, 

^f^mr^ = -^VMr,t) + VUr,t)^ir,t)+gmr,t)\'^ir,t) , 
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with the coupling strength g = Anti^a/M associated to a contact interaction 
with scattering length a for particles of mass M \iAio^l\ lPit61[ [DGPS99J. 

The situation is different when the atoms are trapped by means of an 
electric field interacting with their electric dipole moment, d. Even if alkali 
atoms have ground states with vanishing d, when they are placed in an 
electric field, they acquire an induced electric moment d 7^ 0. Its value is 
proportional to the electric field, with the proportionality constant given by 
the polarizability of the atom. In a first approximation, the polarizability 
is independent of the magnetic quantum number m for atoms of the same 
hyperfine spin / |Ho98j . Therefore, it is possible to store atoms in all the 
states I/; m = — /, ■ ■ ■ , /) in an optical trap. 

The experimental procedure is typically as follows: First, the atoms are 
trapped with a magnetic field. Then, a set of lasers is switched on, defining 
an optical trap superimposed to the magnetic one. Finally, the magnetic 
trap is switched off, and the atoms remain in the optical potential. The 
first realization of such a scheme was done by the group of W. Ketterle at 
MIT in 1997 fMIT98aJ. The remarkable de gree of experimental control on 
lasers allows for a great capacity to modify the trap in a controlled way. In 
particular, it permits the production of very elongated traps, where quasi- 
one-dimensional atomic systems are realized. 

The ability to build such optical traps led the Ketterle group to the 
production of the first spinor condensate with ^^Na atoms [MIT98Jij. The 
ground state of ^^Na has total spin / = 1 and the atoms can be in three 
hyperfine states, m = —1,0,1. In the first stage of the experiment, atoms 
in the \f,m) = |1,— 1) state — which is low magnetic field seeking — are 
magnetically trapped. Then, they are transferred to the optical trap, where 
arbitrary populations of the three states are prepared using radio-frequency 
transitions (Landau-Zener sweeps) |MIT98a.j . Recently, similar experiments 
have been accomplished with ^*^Rb atoms in Hamburg |Hambfl4a] . Georgia 
Tech |(;aTe04| . Gakushuin |(;akun4| and Berkeley |U(]R05| . 

It is interesting to note also that spinor condensates are closely related 
to the pseudo-spin- 1/2 systems realized in a purely magnetic trap by the 
group of Eric Cornell at JILA since 1997 |.IILA97t I.IILA98al I.IILA98hj . In 
these experiments, ^^Rb atoms where trapped in the state \f = l,m = — 1) 
and coupled via a two-photon transition to the state |2,1). Even though 
the states in the / = 2 manifold are metastable, the lifetime of the system 
(a few seconds) is long enough to perform measurements. In such systems, 
remarkable results as spin waves at T > Tc |.TILAfl2] and decoherence ef- 
fects |.TILAfl3h) have been observed, as well as interlaced vortex lattices with 
orthorombic symmetry |.nLAn4c] . This surprising behavior, which contrasts 
with the 'usual' hexagonal Abrikosov lattices encountered in one-component 
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atomic BECs, had been foreseen among others by Mueller and Ho |MHfl2] 
and by Kasamatsu et al. |KTUfl3| . Studies on vortex and spin textures 
for / = 1 spinor condensates can be found in |Muen4j and |RvLSRfl2] . 
Other properties of two-component Bose-Einstein condensates studied in the 
literature are: collective modes |PB98[ ISCn3| IKTUn4j . spin structures in 
the ground and vortex states [HSM IFSM IRDn5| . solitons [Cmn iBAnTl 
lKNF+n4[ IRuo04[ IHern5| . formation of spin domains |KTn4| and Josephson- 
like oscillations |()S99I ILdPT^OOl l(;PS02| . We finally quote also the exper- 
iments on cold gases with two different bosonic atoms by the group of M. 
Inguscio and G. Modugno at LENS, who succesfully created a BEG of "^^K 
atoms by sympathetic cooling with ®^Rb |LENSni] : other mixtures under 

i™n5| . and Li-Gs |Heidni| . 



study are Rb-Gs |Pis 

The theoretical study of integer-/ systems in the context of ultracold 
gases was initiated by T.-L. Ho |Ho98| . Ohmi and Machida |OM98j and 
Law et al. |LPB98| . The correct treatment of the spin degrees of freedom 
of the individual atoms requires to consider the vectorial character of the 
field operator in spin space. In the mean-field approximation, this feature is 
transferred to the order parameter, 



( \ 



(5.1) 



Alternatively, this can be interpreted as having a different order parameter 
for each spin projection. 

An additional peculiarity of these systems is that spin-exchanging inter- 
actions allow for a transference of population between the different compo- 
nents m = — /, ■ ■ ■ , /, subject only to the conservation of the total number 
of particles Af = Nm and of the magnetization of the system, defined 
as Ai = X Nm, being the number of atoms in the m component. 

For instance, two atoms in the state |/ = l,m = 0) may collide and become 
two atoms in the states |1, +1) and |1, —1). Thus, the structure of the exact 
ground state, as well as the dynamics, will depend on the spin-dependent 
part of the Hamiltonian. In fact, it is the importance of the spin degree of 
freedom in these multicomponent systems that has awarded them the name 
spinor condensates, even if the term 'spinor' would not be rigorous from a 
purely mathematical point of view.* 



*An easy introduction to the mathematical concept for 'spinor' can be found in |Wikflfi| . 
More formal texts are |CaT81| for the mathematically-oriented reader, and jCor55llPR,R4[ 
IMF53j for the physics-oriented one. 



74 



Dynamics of spin-1 condensates at finite temperatures 



5.2 Formalism 

As usual for ultracold gases, the low density and temperature of these systems 
allow to substitute the interatomic potential by a contact pseudopotential, 
characterized by its scattering length. To take into account the fact that a 
collision between two spin-/ atoms can happen in any of the channels of total 
spin F = 0, 1, 2, ■ ■ ■ , 2/, this interaction is now written as 

2/ 

V = 5{r^~r2)Y,gFVF, (5.2) 

where Vp projects the two-particle state onto the subspace of total spin F, 
in which the coupling gp is related to the scattering length ap by gp = 
Anh^ap / M . Symmetry constrains F to take even values for identical bosons, 
and odd values for identical fermions. Denoting by fi the spin of atom i = 1,2 
in a collision, and using the identity -F^ = (/i + /2)^, we have 

F=0 

where the summation over F takes into account the possibility of the colli- 
sions occurring through different spin channels. This expression will be useful 
to substitute (some of) the projectors in Eq. ()5.2p by products of atomic spin 
operators. 

We will for definiteness present the results for the case of bosons with spin 
/ = 1 (for instance, ^^Na and *^Rb in their electronic ground state), while 
more detailed calculations for both / = 1 and / = 2 are given in Appendix IbI 
In this case, the previous expressions together with the normalization of the 
projectors, J^f'^f = "Po + ^2 = 1, allow us to write the potential as 

V = Co + C2fl- f2, 

with the coupling strengths given by 

9o + 2^2 47r^2 
Co = ^ = [ao + 202) , (5.3a) 

g2-go 47r;i^ 

C2 = — - — = (02 - ao) . (5.3b) 

Thus, the interaction has a spin-independent part (oc Cq) and a contribution 
(oc C2) that couples atoms with different magnetic quantum numbers m, 
keeping constant the total value M = mi + 7712 for the pair. We remark that 
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this potential is invariant under rotations i?(Q;,/5, 7) in spin space, where 
7} denote the Euler angles defining the rotation (see FigEHl). 
Once the potential is defined, we can write down the Hamiltonian in 
second-quantized form: 

+ |«;„*Jf'»i ■ -f;-*!**} , (5.4) 

where \E'm(^*) (^m(^)) is the field operator that annihilates (creates) an atom 
in the hyperfine state |/ = l,m = 1,0,-1) at point r, and summation 
over repeated indices is to be understood. The trapping potential Kxt('") is 
assumed harmonic and spin independent. Finally, F denotes the vector of 
spin-1 Pauh matrices |Ho98l IZYY04j 

^/010\ 1/° f ^ ^ ^ \ 

= — 1 1 , = -1 1 , = 

\^ 1 / ^v'^ \^ -1 / \ -1 / 

5.2.1 Ground state of the homogeneous system 

Let us first consider the case of a homogeneous system from a mean-field 
point of view: we substitute the field operators by their expectation 
values ipm- Let us also define the spin components of the mean-field order 



Figure 5.1: The Euler angles {a,/5,7} that define a rotation from the 
coordinate system xyz to x'y'z'. 
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parameter through 




(5.5) 



where n is the total density of the system and l^mP = 1- The general 
form to describe a state of a homogeneous spin-1 system with magnetization 
M is 



V 



2 



where 6 and are arbitrary phases, and |,^o| € M is a free parameter. 
The interaction part of the Hamiltonian reads 

Hint 



kl 



Therefore, the spin structure of the ground state is determined by the sign 

of C2 OC 02 — Oq: 



C2 > 0: the energy is minimized by setting {F) = 0, and the system is 
usually called 'antiferromagnetic' |Ho98j . An analogy with condensed- 
matter antiferromagnetism is to be taken with care as there the inter- 
action /i ■ /2 is typically between spins in different positions of space, 
which implies a spatial ordering of them, while here the interaction is 
local. An example of atom with 'antiferromagnetic' interactions is "^^Ua 
in the / = 1 manifold. 

This structure for the ground state is achieved in the state = (0, 1, 0) 
(T denotes transpose) or more generally, due to rotational symmetry 
in spin space, in any of the states 

/ -^e-^°sin/3 
UF = e'n cos (3 

\ ^e^" sin (3 

where 6' is a global phase that does not affect any expectation value, 
and {a,/?} are the Euler angles that fix the quantization axis. This 
ground state presents a U{1) x S"^ symmetry and is usually referred to 
polar ground state |Ho98j . 
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C2 < 0: the energy is minimized by setting \ {F) \ = 1, and the system 
is called 'ferromagnetic' [Ho98|. The ground state configuration of the 
spin components reads 




e"^" cos^ f 
\/2 cos I sin I 
e^" sin^ I 



Note the coupling of the Euler angle 7 to the global phase 9. The 
symmetry is now SO (3) |Ho98| . 

A similar treatment follows for larger spin values, the main difference 
being the greater number of scattering parameters that enter into the play 
and, thus, the richer phase diagram that emerges. For instance, the case 
of / = 2 has been studied by Ciobanu et al, who show that there are 3 
possible phases for the ground state depending on the values of the couplings 
Co oc 4a2 + 804, Ci oc 7ao — I0a2 + 804 and C2 oc 04 — 02 ICYHOOj (see 
also (KMI)- 

It is also worth noticing that the presence of a homogeneous magnetic 
field can drastically affect the structure of the ground state, as was analyzed 
by Zhang et al. for the case of / = 1 |ZYYn3| and for / = 2 by Ciobanu 
et al. |(^YHnn j and more recently by Saito and Ueda [SUOBb] . Also the 
dynamics of the system is strongly infiuenced by external magnetic fields as 
studied theoretically for the case of spin 2 by Saito and Ueda |SUn5aj . and 
shown experimentally by the Hamburg |Hambfl4b] and Gakushuin [Gakun4j 
groups. This notwithstanding, in order to reduce the number of parameters 
entering our simulations and thus have a clearer picture of the effects of 
temperature, we will restrict our analysis to the case of a vanishing magnetic 
field. 



5.2.2 Dynamical equations and transfer of population 

We focus now on the dynamical evolution of a trapped / = 1 condensate 
and, in particular, on the possible influence on it of thermal effects. In 
experiments, practically any initial configuration N{t = 0) = {Ni, Nq, N_i) 
can be produced,^ and it is interesting to analyze e. g. whether the system 
is expected to converge in the time scale of an experiment to its ground- 
state configuration, or the spin structures that can appear in the dynamical 
process. 

'^Here AT is to be understood just as a shorthand for {Ni,No, N^i), and not as a true 
vector in spin space. 
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The equations of motion for the mean-field order parameters can be de- 
rived from the energy expression in a mean-field formulation, 

£[^] = J dV (-^V^ + Kxt) i^m 

+ ^rmrji^ji^m + ^rmr^F^k ■ F^iMk^ , (5.6) 

by functional differentiation according to ihdilim/dt = SS/Sip^. One readily 
obtains (see Appendix m for details) |PT:R+99[ IZYYOH] 



in- 



2M " 



dt 

where we have defined the transference terms 



i^m + C2T^ , (5.7) 



(5.8) 



which render these equations different from the usual number-conserving 
Gross-Pitaevskii equation for one-component condensates. In fact, the C2- 
term in the Hamiltonian does not couple directly the different \f,m) states 
in which one atom can be found — as the radio-frequency field couples the 
|1, —1) and |2, 1) states in the JILA experiments — , but the two-atom states 
\Z) := |/i,mi)|/2,m2) = |1,0)|1,0) and \U) := |1,1)|1,-1) (properly sym- 
metryzed). The effective potentials that will determine the spatial dynamics 
of each component read 

V±f =Kxt + Con + C2 (±^1 + no =F n_i) , . 

K)'^=Kxt + Con + C2(ni + n_i) ' ^^"''^ 

with rimir) = \iprn{f)\^ being the density of atoms in the |l,m) state and 
n(r) = J2m ^m{f) the total density, normalized to the total number of atoms 
M. The population of the hyperfine state |l,m) is = f d^r iV'mP, and 
Ni + Nq + = M is a constant of the motion. 

In analogy to what happens for the Gross-Pitaevskii equation of a spin- 
polarized condensate IDGPSQO] . equations (I5.7jl can be rewritten in the form 
of continuity equations, but now a balance term accounting for the transfer 
of populations between the components will appear on the right-hand side. 
Indeed, with the usual definition of the quantum current applied to each 
component, 

h 
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one readily finds 

-^ + '^ ■jm = ShUr,t), (5.10) 

where 6hm{r,t) = —{2c2/h) Im[Tm?/'m] is the transfer of populations between 
spin components per unit time. 

The dynamical equations for m = ±1 are invariant under time-reversal 
symmetry; therefore, 5hi = 6h_i, and the conservation of magnetization 

implies Shq = -26n±i = {2c2/ih) {ipiipl^ip^i - ^/'J'^'oV'-i)- 



5.3 Numerical procedure 

We have studied a system of ^^Rb atoms in their electronic ground state 
manifold / = 1, trapped in a very elongated trap as in the Hamburg experi- 
ment |Hambn4b] . In particular, we have considered 20 000 atoms confined in 
a trap characterized by frequencies u± = 891 Hz and = 21 Hz, which re- 
sults in an aspect ratio uJz/uj± = 0.024 <^ 1. For this atom, 02 = 100. and 
ao = 101. Sob |vKKHV02] . which results in C2 < 0, and the expected behav- 
ior is 'ferromagnetic'. Moreover, we have simulated a purely one-dimensional 
system in a trap of frequency Uz, as it has been shown that elongated spin-1 
systems can be safely considered one-dimensional for most purposes |ZY05j . 

We have solved Eqs. (I5.7jl and analyzed them in terms of the transference 
terms (I5.8jl and the effective potentials (j5.9jl . To this end, at each time step 
dt we have split the evolution in two parts: the kinetic energy part, and the 
rest of the Hamiltonian. The evolution due to the kinetic energy has been 
solved by transforming the wave function into momentum space by a Fast 
Fourier Transform subroutine, and then the operator exp[—ip'^/{2m)dt] has 
been applied to it. Finally, the wave function is transformed back into real 
space, where the effect of the trapping potential and the coupling terms has 
been taken into account by means of a fourth order Runge-Kutta algorithm. 
This splitting of the full evolution in two parts introduces, in principle, higher 
order terms when the separate parts do not commute with each other [see 
Eq. (I6.12|) ]. However, for short enough time steps dt, these terms can be 
neglected. To check the accuracy of our procedure, we also solved the evolu- 
tion entirely with the Runge-Kutta algorithm, computing the gradient terms 
numerically in r-space: both methods give the same results for the test cases 
analyzed, but the split method can make use of much larger time steps and, 
thus, calculations are faster. 

We consider that initially a quasi-pure condensate in the m = spin com- 
ponent is populated. Numerically, spin mixing requires at least a small seed of 
atoms populating the other components, but we keep the total magnetization 
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equal to zero for the following reason: It has been shown that, for the case of 
zero magnetization, not all the phases 6m in Eq. affect the spin dynam- 
ics, but it is governed by only one relative phase Q = 260-61-6^1 |PT:R+99j . 
Therefore, we set 6*1 = 6'_i = and vary only 6q. This relative phase is a very 
important parameter, as it can freeze completely the spin dynamics, or make 
it faster or slower [PLR+QQ] . In the experimental scheme, the phases are not 
well controlled, and vary from shot to shot. Thus, in order, to reproduce 
experimental results, we will randomly draw 20 values for 6 in the range 
[0,27r), and make an average of the corresponding results. 

5.4 Dynamical evolution at zero temperature 

First of all, we perform several simulations at zero temperature, to be com- 
pared with the work of Yi et al. |YMSYn2] . These authors work in the 
so-called Single Mode Approximation (sma) which, based on the fact that 
|c2| -C Co, assumes imif) = Cm [PLR+QQ] . This approximation allows for 
some analytical insight into the properties of condensates with spin degree of 
freedom, and has therefore been extensively used. However, its accuracy was 
analyzed by Pu et al. |PLR,"'"99] . who established that it is only adequate for 
small particle numbers, when the contribution from the C2 term is negligible. 
As M increases, the contribution from the spin-exchange term is more and 
more relevant and the SMA becomes invalid. In our system, Af = 2 x 10'^ 
and we have checked that indeed the density profiles for the different com- 
ponents cannot be transformed into one another by a simple rescaling, i. 
e. nm{r)/nm'{r) ^ constant Vr (m 7^ m') (see Figs. 15.31 and 15.61 below). 
Therefore, we will not use the SMA in our work. 

The results of our simulations at zero temperature are summarized in 
Figs. I5.2H5.4I In figure 15.21 we plot the population of each spin component 
as a function of time, for the initial populations {Ni/M, Nq/M, N_i/J\f) = 
(0.5%, 99%, 0.5%). Dashed lines correspond to the dynamical evolution with 
an initial relative phase 9 = 0, and solid lines stand for the average over 20 
random initial relative phases. In absence of a magnetic field gradient the 
dynamical evolution of the m = ±1 spin components is the same, and the 
curves Ni{t) and iV-i(t) coincide; besides J2^mit) = M is always fulfilled. 

For the case of fixed 9 = (dashed lines) , we observe several oscillations 
of the populations as a function of time, with a period r ~ 50 ms. These oscil- 
lations are reminiscent of Josephson oscillations between the two-atom states 
|1, 0)1 1,0) and |1,1)|1,— 1), coupled by the C2-term of the Hamiltonian, but 
appear to be damped. We have checked that these oscillations are coherent 
in a homogeneous system, i. e., after a certain period tj ^ h/{c2n) ~ 50 ms 
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(n = 10^^ one has N{t = tj) = N{t = 0), while at t = rj/2 the popu- 

lations are reversed, N±i{t = tj/2) = No{0)/2 and No{t = tj/2) = 2iV±i(0). 
In the case of the trapped gas, the inherent nonlinear character of the in- 
teracting system, together with the discreteness of its spectrum, induce a 
'collapse' |ENS8ni IM(]WW97| or 'dephasing' |VT:99[ [LegOlt 101:021 of the 
oscillations. Moreover, as the period of oscillation depends strongly on the 
relative phase 0, an average over various values induces a further reduction of 
the amplitude of the oscillations, as shown by the solid lines in the figure. At 
much longer times, it has been shown that the full quantum solution presents 
revivals of the exchange of populations |PLR"'"99| . Nevertheless, this result 
cannot be retrieved in a mean-field calculation such as ours IMCWWW) . 

In both our calculations (with or without average), the magnetization 
A4 = J2m ^ IS couscrvcd duriug the time evolution, as it has been ex- 
perimentally observed |(]aTen4| . The oscillations between the populations 



Figure 5.2: Population of the spin components as a function of time for 
the initial configuration {N^/M ,Nq/M ,N_i/ M) = (0.5%, 99%, 0.5%) and 
T = 0: the light, green lines stand for Noit), while the dark, black lines 
are N±i(t). The results obtained with one initial phase 6 = are denoted 
by dashed lines, while solid lines show numerical results averaged over 20 
random initial relative phases 0. 
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of the m = and m = ±1 states are not fully coherent but present a dy- 
namical instability around t ~ 100 ms, when the large amplitude oscillations 
become small amplitude oscillations |SU05b| . At this point, the popula- 
tions start to oscillate around the ground state configuration of the system 
with zero magnetization which, in absence of an applied magnetic field, is 
(25%, 50%, 25%) |ZYY03| . We note that these numerical results are in qual- 
itative agreement with the experimental measurements of Ref. |GaTe04j ob- 
tained in a strongly anisotropic disk-shaped trap, where the relaxation to the 
steady state is also not monotonic but presents a few damped oscillations. 

To further understand the spin dynamics, we now analyze the time evo- 
lution of the density profiles of the diff'erent spin components in the trap, 
in search for the possible existence of spin waves or the formation of spin 
domains. We plot the evolution corresponding to the run with fixed 9 = 
in Fig. 15.31 At the initial stages of the evolution, the population of the m = 
component decreases due to the spin-exchange interaction, and the ±1 spin 
components start to be populated by equal amounts, keeping the magnetiza- 
tion of the initial state. In fact, the magnetization is locally conserved along 
the time evolution, a fact not required by the symmetries of the Hamiltonian. 



Figure 5.3: Density profiles of the spin components at the times indicated 
(in ms) and T = 0, with the same color labelling as in Fig. 15.21 The initial 
configuration corresponds to N{t = 0)/Af = (0.5%, 99%, 0.5%), and 6 = 
(dashed lines in Fig. 15. 2j) . 
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For times t < 100 ms, the conversion of atoms from m = to m = ibl 
occurs mainly in the central part of the condensate, where the density is 
higher and thus the coupling among the different spin components is more 
effective, see Eq. ()5.9p . Later, at the time when the dynamical instability sets 
in, the ±1 spin components swing back to the component, giving rise to the 
formation of a spin structure with small domains. The m = and the m = 
±1 domains are miscible, forming what can be named as 'mixed domains'. 
This fact is related to the 'ferromagnetic' character of the interaction: if 
in the m = ±1 domains the spin is locally oriented (anti) parallel to the 
quantization axis, and in the m = domains it is perpendicular to it, in 
these 'mixed domains' the spin is locally oriented along a direction neither 
fully parallel nor perpendicular to the quantization axis. Thus, the succession 
of domains should reflect a smooth twist of the spin along the system, as also 
suggested by |SU05h|. 

The identical location of the m = ±1 domains is due to the symmetry in 
the corresponding initial proflles, together with the fact that the dynamical 
equations are symmetric under the relabelling 1^—1. Note also that the 
number of small spin domains does not grow indeflnitely, but it is limited 
by a characteristic size /dom that depends on the internal coupling between 
different spin components. We remark that the total density proflle is con- 
stant during all the simulation, as expected for a trapped condensate in the 
Thomas-Fermi regime. Indeed, in this regime the spin-independent inter- 
action ~ cqu ~ 100 nK is dominant over the kinetic energy terms, thus 
preventing the formation of total-density modulations |HDR"'"Ol| . 

It has been argued that the appearance of the small domains is intimately 
related to a dynamical instability of the system. Therefore, the time scale of 
this process can be estimated from an analysis of the spectrum of the sys- 
tem. For two-component condensates, Kasamatsu and Tsubota |KT04| have 
observed a very fast appearance of structure in the spin density distribution 
of a ^^Na system. They estimate the time for the growth of domains as 
T"growth = 27r/|f2_|, where f2_ is the excitation frequency with largest imagi- 
nary part. In their case, Tgrowth ~ 25 ms. 

A similar analysis for / = 1 condensates has been recently presented 
by Zhang et al. |GaTe05bj . For the case of 'ferromagnetic' interactions, the 
instability is expected to emerge at times ~ 2'7ih/{\c2\n). For ^^Rb at densities 
n ~ n{z = 0) ~ 5 X 10^^ cm~^, this estimation is around 80 ms, which is in 
fair agreement with our results, cf. Fig. 15.31 They also estimate the typical 
size of a spin domain at long times as Adom = '^T^/hnst, where fcinst is the 
highest momentum for which an unstable mode exists [i. e., lmu{ki^st) 7^ 0]. 
With this formula, they get Adom ~ 13 /im, a value consistent with our 
results. For 'antiferromagnetic' systems such as / = 1 ^^Na, they predict no 
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instability and, therefore, no domain formation, a point that we will address 
in Sect. 15 .5 .31 

It is possible to interpret the dynamical evolution of the various spin 
components in terms of the effective potential V^{z,t) felt by each one and 
through the continuity equations To this end, we plot in Fig. 15.41 the 

effective potentials V^^^ and = Kff [top panels], and the local transfer 
of population dtUQ^z^t) and dtni{z,t) = dtn_i{z,t) [bottom panels] at times 
t =0, 40 and 80 ms, for the same initial conditions as in Fig. 15.31 As is 
clear in Eq. ()5.10|) . dtUm has two contributions. For a fixed time, 6hm{z,t) 
represents the local variation in the occupation of the |l,m) state due to 
spin-exchanging collisions, while —V ■ jm takes into account the change in 
rim due only to motion of m-atoms. At the initial stages of the evolution, 
t < 40 ms, the densities of the m = ±1 components are very small and the 
dominant contribution is for all components. For example, at t = 40 ms 
6ho{z,t) < and consequently 6hi{z,t) = 6h_i{z,t) > 0: the C2-term in the 
Hamiltonian is favoring the conversion of |1, 0)|1, 0) pairs into |1, 1)|1, —1), as 



Figure 5.4: Effective potentials V^{z,t) in kHz [top] and population 
transfers 6hm{z,t) in /im~^-ms~^ [bottom] at T = for the initial configura- 
tion N{t = 0)/Af = (0.5%, 99%, 0.5%) and 6 = (dashed lines in Fig. EH, 
at times t = 0, 40 and 80 ms. 
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seen in Fig. 15 .21 for the first stages of the evolution. At t = this is also true, 
but it cannot be appreciated from the scale of the figure. Moreover, since 
the minimum and maximum of 6nQ{z,t) and 6h±i{z,t), respectively are at 
the center of the condensate, the spin-exchange mainly occurs at the central 
region, as we have already observed in Fig. 15.31 

At t = 80 ms Sn±i is positive at the boundaries of the condensate, and 
negative in the central region, whereas Suq has the opposite behavior. There- 
fore, the population with m = ±1 is converting into m = atoms at the 
center while still increasing at the boundaries. In fact, at this time, the total 
tranfers ANm{t) = J dz 6hm{z,t) between m = and m = ±1 almost bal- 
ance, which is reflected in the vanishing slope of the curves Nm(t) of Fig. 15.21 
At this stage of the evolution, the population of the m = ±1 states starts 
to be relevant, and the contribution of the currents is no longer negligible. 
Typically, it has the same sign as the spin-exchange contribution, but it is 
smaller in magnitude. 

These features can also be understood in energy terms from the panels 
showing the effective potentials. For example, at t = the higher energy for 
the m = atoms induces them to convert into m = ±1 pairs, especially in 
the center of the cloud, where the difference V^^''^ — V^^ is larger. At time 
t = 40 ms, V^''^ = at the center and, therefore, the transfer of atoms 
will slow down there and reverse its sense shortly after, while the population 
of the ±1 components will continue to grow on the sides of the cloud (cf. 
Fig. I5.3|l . Later on, at t = 80 ms, the energy of the m = component is still 
larger on the boundaries than that of the m = ±1 components, but it has 
become smaller in the center of the trap; therefore, 0-atoms will continue to 
convert into ±l-atoms in the boundaries, but the reverse will happen at the 
center, as confirms the lower panel. 

5.5 Dynamical evolution at finite temperature 

5.5.1 Introducing temperature fluctuations 

We would like now to analyze if thermal effects may have some noticeable 
effect on the dynamics of spinor condensates, making their expected evolution 
differ from the results found in the previous section for T = 0. 

For single-component condensates, temperature effects on the conden- 
sate fraction, dynamics, and damping of excitation modes have been studied 
extensively. Usually, at low enough temperatures, thermal excitations can 
be accounted for within the Bogoliubov-de Gennes (Bdo) |dG66[ INP90| or 
Hartree-Fock-Popov (hfp) |ZNG99| frameworks. Recently, this last scheme 
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has been applied to study finite temperature effects in the equilibrium density 
distribution of the condensed and non-condensed components of / = 1 atoms 
optically trapped |ZYYn4| . In our work, we will make use of the BdG theory 
to describe the thermal clouds present in multicomponent condensates. 

For a highly elongated one-component system, it has been shown that 
at low temperatures thermal fiuctuations are relevant in the phases of the 
field operator, while they can be disregarded on the density profile of the 
ground state |KK65I IR,C65j . To be precise, three regimes can be iden- 
tified [PSWOnj . according to the value of the temperature relative to the 
degeneracy temperature T^eg = Mhoj |KvD9 61 and a critical temperature for 
phase fiuctuations Tg = Tdeghw/fi- (a) a true condensate (density and phase 
fiuctuations are small in the ground state) at T -C T^; (b) a quasiconden- 
sate for Tg <^ T <^ ^dcg: the density has the same profile as for the true 
condensate, but the phase fiuctuates on scales smaller than the cloud size, 
and the coherence properties of the phase are drastically modified; (c) finally, 
for T > Tdeg, both phase and density fiuctuate, and the system is no longer 
quantum degenerate. 

In the case of a one-dimensional, harmonically trapped single-component 
condensate, the Bogoliubov-de Gennes equations for the low-lying excita- 
tions, 



^i^i = ( ^^^^ + ^^^t - l^ + '^an] Uj + griQVj (5.11) 



-— + Kxt - /i + '2^gn ) Vj + gnoUj 



(5.12) 



can be solved exactly |Str98] . obtaining the corresponding energy eigenvalues 
tj = ^/ j{j + l)/2huj and wave functions. These can be written in terms of 
Legendre polynomials as^ 



2j + 1 /i 
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TF 



1/2 
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TF 



with the chemical potential given by /i = ;k^(37Va/V32)2/3, where a gives 
a measure of the strength of the interactions as compared to the trapping 
potential energy |PSWnfl| . and i?TF = \j2^/m/uj is the Thomas- Fermi radius 
of the initial condensate. 



■I-For quasi-lD condensates, the phase fluctuations at T = are described by Jacobi 
polynomials in contrast to exactly ID condensates (see jPSWflfll IStr98j l. Nevertheless, 
already at T ~ 0.2Tc quantitative difference between thermal e ffects in quasi-lD and 
strictly ID cases becomes irrelevant (compare |HDR + Oll iKSS+na] !. 
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As the density fluctuations with respect to the equilibrium profile Ueqiz) 
may be disregarded, one can write ^{z^O) = ^ynec^{z) exp[i6{z)] for the par- 
ticle anihilation operator, with the phase operator given by [PSWOO 



A/4neq(^) 



]\z)aj + H.c. 



where hj anihilates a quasiparticle in mode j. Thus, in the temperature range 
To <^T <^ Tdeg, where only phase fluctuations are relevant, one can simulate 
the thermal fluctuations in a one-dimensional condensate by generating a 
thermal phase 9^^{z) according to 
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and adding it to the mean-fleld order parameter before the real-time simu- 
lation starts |Ha,mbni| IKSS+AB] . An analogous treatment for the case of a 
multicomponent condensate shows that, assuming that initially almost all the 
atoms are in a single component, the BdG equations for the non-condensate 
atoms depend mainly on the excited fraction in this most populated hyperflne 
state. Therefore, one can simulate again temperature effect by introducing 
one such fluctuating phase for each spin component, 
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(5.13) 



that will be added to the corresponding 9m of Eq. (jB.ljl . 

Here a™ (a™*) are complex amplitudes that replace the quasi-particle an- 
nihilation (creation) operators in the mean-field approach. In the numerical 
calculation, in order to reproduce the quantum-statistical properties of the 
phase fiuctuations, and a™* are sampled as random numbers with zero 
mean and (la^p) = [exp^ej/ksT) — the occupation number for the 
quasi-particle mode j |Hambni] . Let us note that this stochastic procedure 
will translate in different initial conditions for m = 1 than m = — 1, and will 
therefore break the symmetry between them. 



5.5.2 Results for a 'ferromagnetic' system: ^'Rb 

We have performed a series of simulations following the scheme above to 
introduce temperature fiuctuations in the phases. The dynamical evolution 
of the population of the spin components that we obtain is plotted in Fig. 15.51 
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Figure 5.5: Population of the spin components as a function of time for 
the initial configuration N{t = Q)/M = (0.5%, 99%, 0.5%) and T = Q.2T^. 
The labelling of the lines is as in Fig. 15.21 
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for the case of T = 0.2Tc and the same initial populations as in Fig. 15. 2[ i. 
e., N/M = (0.5%, 99%, 0.5%). Here, = ^^nLU J ln{2^^) is the critical 
temperature of Bose-Einstein condensation for a single-component ID Bose 
gas in a harmonic trap of frequency |KvD96| . As in previous figures, solid 
lines correspond to the numerical results averaged over 20 random values for 
the phase 9, and dashed lines to a single run with 9 = 0. One observes 
that the interaction of the condensate atoms with the thermal clouds smears 
out the oscillations present at zero temperature (cf. Fig. I5.2|l and leads to 
an asymptotic configuration with all components equally populated. We 
note that this spin distribution has also been experimentally obtained from 
some initial preparations with zero total spin |Hamh04h| . This eff'ect can be 
understood in terms of the free energy of the system, T = E — TS\ as the 
spin-related contribution to the energy, ~ ~ 1 nK, is relatively small 
compared to typical experimental temperatures Tcxp ~ 100 nK |MIT98bj . 
the entropy contribution to T dominates, and the equilibrium configuration 
will be the one that maximizes iS, i. e., with equipartition. 

In a spin-polarized, elongated condensate, the occurrence of phase fiuctu- 
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ations due to thermal excitations induces modulations of the density during 
an expansion of the system in a time scale of the order of 10 ms ^SS"'"03J. 
These density fluctuations are not expected inside the trap as the mean field 
prevents them, but in a time-of-fiight experiment they can be used to quan- 
tify the phase fluctuations of the trapped system and also as a thermometry 
probe iKSS+OBj . In a multicomponent condensate, one can follow the same 
argument to show that modulations of the total density should not be ex- 
pected in the trap. 

However, the situation is different for the densities of the various spin 
components. Indeed, the spin-exchange interaction (~ C2n) is much weaker 
than the scalar mean field (~ cqu), and it is not strong enough to prevent 
spin-density fluctuations, which lead to the formation of spin domains, cf. 
Fig. 15.61 This process is much faster at finite temperature than at T = 0; 
for example, the number of spin domains at time t = 80 ms for a temper- 
ature T = 0.2Tc is larger than at time t = 240 ms for T = 0. Such a 
fast domain formation has also been observed in simulations of pseudo-spin- 
1/2 systems [KT04ij . This result contrasts with the relatively long times for 
the emergence of domains predicted by the theory and simulations of Zhang 
et al. |GaTe05b] and obtained in the simulations by Saito and Ueda |SU05bj . 
We think that this discrepancy may be due to the presence of the thermal 
clouds. Indeed, it has been argued that spin oscillations in the thermal cloud 
can considerably affect the spin dynamics in a finite-temperature conden- 
sate |.nLA02| (see also |()L02| ). With regards to this point, it is interesting 
to note that, even though our way of introducing thermal phase fiuctua- 
tions does not take into account possible spin correlations among the clouds, 
the fact that the initial system is almost completely in one component may 
prevent such correlations from being relevant for the later dynamics. 

Moreover, the density profiles for the m = ±1 components at finite tem- 
perature are no longer equal as was the case at T = 0. This is due to the 
existence of different thermal clouds for each component |.TIL AOSbj IHamb04cl 
IHamb04'a] . which is accounted for in the simulation by having different ran- 
dom profiles for the corresponding thermal phases, Of'i^z) ^ Q^\{z). As a 
consequence, the local magnetization ni(z,t) — n_i(2;,t) is no longer con- 
stant as at T = 0. Nevertheless, as expected, the total magnetization is still 
a conserved quantity along the time evolution. 

The typical size of the domains /dom can be estimated from our data to 
be ~10 /xm. Ueda |Ued01| and Saito and Ueda |SU05b| have studied the 
excitation modes of the system within a Bogoliubov-de Gennes scheme, and 
determined that it should undergo a dynamical instability through modes 
related to spin waves which carry angular momentum The most unsta- 
ble modes according to them are those with a momentum fcinst that satisfies 



90 



Dynamics of spin-1 condensates at finite temperatures 



Figure 5.6: Density profiles of the spin components at different times (in 
ms) at T = 0.2Tc for the same initial configuration as before (dashed line 
in Fig l5.2|) . Again, the total density keeps the same profile throughout the 
simulation. 
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2vr/A;inst ~ ^dom- Their simulations (for ID and 2D systems) also present this 
typical domain size |SUfl5h| . They also note that this size should depend on 
the presence of a homogeneous magnetic field: for larger fields, the domain 
size would grow, reaching ~ 16 /xm for -B = 1 G. This value is closer to the 
experimental results |GaTen5a] than ours. Also a recent work in an elon- 
gated 3D system by Zhang et al. |GaTefl5h] reports a value similar to the 
experimental one including in the simulations a magnetic field B = 0.3 G. 
Nevertheless, it is also possible that other thermal effects might be impor- 
tant. For example, spin excitation modes should be taken into account in 
our simulation on an equal footing to the phase fluctuations. To this end, 
the wave functions corresponding to the spectrum of spin excitations found 
in |Ho98| ISUn5b| IGaTeOBb] need to be calculated. Then, they should be 
thermally populated at the beginning of the simulation in a way similar to 
what we have done for the phase fluctuations. 

Finally, we note that the results obtained by Saito and Ueda are for 
a system at zero temperature. The reason why they obtain a dynamical 
separation of the ±1 domains — contrary to our results, see Fig. 15.31 — is 
that their starting input profiles have a small magnetization, which triggers 
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the formation and development of spin waves along the system. Also, the 
dynamical instability appears in their calculations at later times t ~ 300 ms, 
even though they have more particles in the system, M = 4 x 10^. These 
differences are most probably due to the fact that the inhomogeneity in the 
initial conditions that we impose through the phase fluctuations is more 
important, and favors the appearance of domains. 

5.5.3 An 'antiferromagnetic' system: ^'^RbAF 

After studying the system of '^''Rb (/ = 1), that is 'ferromagnetic', we turn 
our attention to an 'antiferromagnetic' case, to check the prediction of Zhang 
et al. that for such systems no dynamical instabilities occur and, therefore, 
no domain formation is to be expected |GaTe05b] . A system naturally 'an- 
tiferromagnetic' is ^^Na in the / = 1 hyperflne manifold. However, in order 
to do a direct comparison with the results of the previous section, we prefer 
to work in a model system, which is formed by ^^Rb atoms with / = 1 but 
where we have changed the sign of the coupling constant C2. In this way, 
only the spin dynamics is expected to change, while nothing different should 
happen to the total density as compared to the previous results. We will 
denote this artiflcial version of rubidium as ®^RbAF- 

As before, we have performed simulations at zero and finite temperature. 
The results are presented in Figs. 15. 7H5. 81 and they are readily seen to be quite 
different from those corresponding to the ferromagnetic case. In figure 15.71 
we present the evolution of the population of the different spin components 
at zero and finite temperature. The zero-temperature results, displayed by 
dashed lines, consist in almost perfect oscillations between the initial pop- 
ulation N/N = (0.5%, 99%, 0.5%) and the exact mean-field ground state 
(50%, 0%, 50%) |Z YY03j . This behavior can be interpreted as Rabi-like os- 
cillations of pairs of atoms between the states \U) and \Z), defined above [see 
discussion after Eq. (j5.8jl ]. Indeed, the ground-state configuration is basically 
formed by J\f/2 pairs in state \U) and the initial one has the same number of 
\Z) pairs; moreover, both states satisfy \{F) \ = and are almost degenerate. 
However, the exchange of populations is not complete, and the overall data 
cannot be fitted by a simple sinusoidal function, but the oscillations seem 
to 'speed up' as time goes by: the two first maxima of Nq/M are separated 
~ 230 ms, while the second and the third are ~ 215 ms apart. Again, this 
fact should be attributed to the nonlinearity of the sytem. 

The results for the simulation at T = 0.2Tc (solid lines) present a very 
different behavior. For the case = 0, there is only a very slow decay of the 
population of the m = component, with the corresponding increase of the 
number of atoms in m = ±1. This is due in part to the fact that the chosen 
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Figure 5.7: Population of the spin components as a function of time for 
the initial configuration N/Af = (0.5%, 99%, 0.5%) and 9 = for the 'anti- 
ferromagnetic' system ^''RbAF (see text). The dashed lines are the results at 
T = while the solid lines have been calculated at T = 0.2Tc. 




Time [ms] 

starting configuration is energetically very close to the ground state for the 
antiferromagnetic system, which will result in a long thermalization time. We 
also remark that in neither the zero-temperature or the finite-temperature 
steady-state is observed in the simulated time of 500 ms, contrary to 
what happened for the 'ferromagnetic' case (cf. Figs. 15.21 and 15 .5^ . One may 
expect, however, to converge to an equilibrium state with longer simulations. 

The very slow evolution of the spin transfer in the finite-temperature case 
is also refiected in the density profiles, that we present in Fig. 15.81 Again, 
solid lines stand for the T = 0.2rc case, while dashed lines are for the T = 
calculation, and the color coding is the same as in previous figures. At zero 
temperature, in the initial stages t < 80 ms, we see an evolution resembling 
that of the 'ferromagnetic' ®^Rb (compare Fig. l5.3|) . However, for t = 240 ms 
the system practically turns back to its initial distribution, a behavior very 
diff'erent from the 'ferromagnetic' one, where one could observe the formation 
of spin domains. The diff'erences are still more pronounced for the simulation 
at finite temperature, where there is essentially no evolution beyond some 
small fluctuations, coming from the random thermal phases. Only at very 
long times t ~ 500 ms (not shown in the figure for clarity) the m = ±1 clouds 
start to be visible on the scale of the figure. 
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Figure 5.8: Density profiles for the 'antiferromagnetic' system ^'RbAP for 
times t = 0, 80, 160, 240 and 500 ms. The dashed lines are the results at 
T = while the solid lines have been calculated at T = 0.2Tc. The initial 
configuration is N/M = (0.5%, 99%, 0.5%) with 6 = 0. The m = profile 
is plotted in green and the m = 1 profile in black, while the red line in the 
last plot is for the m = —1 component. 




Position Position [j^tn] 



Compared to the 'ferromagnetic' case in Fig. 15. 6[ one remarks the dif- 
ficulty in forming separate spin domains, even though the thermal phases 
are different for all components. This must be attributed to the 'antiferro- 
magnetic' character of the system under study, that tends to keep atoms in 
states m = ±1 together to minimize the local spin, \{F)\{z) = 0. On the 
other hand, one expects a separation of m = ±1 atoms from m = atoms, as 
observed experimentally for the case of ^^Na |MIT98b| . which also explains 
the difficulty in populating the m = ±1 components starting from an almost 
pure 0-condensate. 



5.6 Conclusions 

In this chapter, we have studied the spin dynamics of a spin-1 ^^Rb conden- 
sate in a highly elongated trap, which we have modeled as a one-dimensional 
system. We have solved the three coupled dynamical equations for the spin 
components within a mean-field framework without any further approxima- 
tions. This is in fact necessary, since approximated treatments, such as 
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SMA, frequently mask some of the aspects of the dynamics. We have also 
considered finite temperature effects. For one-dimensional systems at very 
low temperatures, these effects manifest themselves basically as phase fiuc- 
tuations, which we have incorporated using the approach of Ref. [HDR+nij 
based on a thermo-statistical population of the excitations of the system. 

We have found that the spin dynamics towards the steady state is not 
monotonic but rather slowly damped, involving a coherent exchange of pop- 
ulation between the various spin components. At finite temperature, the 
coherent oscillations of populations among the m = and m = ±1 states 
are strongly damped. The internal coupling of the spin components together 
with the presence of a different thermal cloud for each component lead to the 
formation of numerous spin domains within the condensate. The size of these 
domains defines a characteristic length scale in the system, /dom ~ 10 /im, 
which does not decrease with time, and seems to be intrinsically determined 
by a dynamical instability of the excitation modes of the 'ferromagnetic' sys- 
tem, without noticeable thermal effects. Indeed, similar simulations for a 
system in which we have artificially changed the sign of the spin-exchange 
interaction present no formation of domains even at finite temperatures, as 
expected from theoretical considerations |GaTefl5b) . 

For a 'ferromagnetic' condensate at zero temperature and with initially 
almost all atoms in the m = component and zero magnetization, the pop- 
ulations of the spin components oscillate around the configuration which 
corresponds to the ground state, (25%, 50%, 25%). On the other hand, at 
finite temperature the interaction of the condensate atoms with their ther- 
mal clouds leads to equipartition in populations (1/3,1/3,1/3). However, 
recent simulations at zero temperature with a very small initial magnetiza- 
tion M-jM = 10""^ also drive the system to stationary states far from the 
ground state. Moreover, in these simulations formation of multiple spin do- 
mains is also observed. Therefore, the true role of thermal effects on the 
dynamics (beyond that of breaking the symmetry of the initial conditions of 
the m = ±1 components) is not clear yet, and deserves further consideration. 

The results for an 'antiferromagnetic' case are very different. In partic- 
ular, the dynamics of populations and density profiles is very slow, and no 
domains appear for times t < 500 ms starting from an almost-pure m = 
condensate. We have explained this result by the 'antiferromagnetic' charac- 
ter of the interactions in this system together with the quasi-degeneracy of 
the initial configuration with the mean-field ground state. 

Our results might also be relevant to the question of decoherence. Our 
simulations, and in particular the finding of very fast domain formation, 
suggest that decoherence is enhanced as the number of components in the 
system is enlarged. Of course, there are many open questions connected to 
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this, e. g., whether the formation of domains goes along with a loss of phase 
relations, and gives rise to some enhanced (generalized) phase fluctuations. 
It seems also interesting to study the mechanism underlying the speeding of 
domain formation by temperature. To this end, one might consider incor- 
porating spin-density fluctuations already at the start of the simulation in a 
way similar to what has been done for the phase fluctuations. 
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Book II 

Two-dimensional helium systems 



Motivation: helium as 'the' 
quantum hquid 



In this second part of the thesis, we aim at studying a variety of hehum 
systems, which are quite different from the ultracold, dilute gases. Hehum has 
for a long time been considered 'the' quantum fluid, as many of its physical 
properties can be directly related to the quantum laws: from the fact that 
helium remains liquid down to zero temperature at saturation vapor pressure, 
to the superfluidity present in both the bosonic ^He and the fermionic ^He 
species. In fact, since the pioneering experiments of H. Kamerling Onnes 
to cool down helium and other substances |Kamll[ IKa,m67| . until the more 
recent claim for the observation of a supersolid phase in helium [KC04J, the 
study of helium has received a lot of attention. This attention has rewarded 
us with a more profound understanding not only of this atomic species, but 
of the implications of quantum laws when applied to a many-body system. 

Among the most striking experiments, one can mention: the liquefaction 
of helium by H. Kamerling Onnes (1908); the discovery of '^He superfluid- 
ity by Pyotr Kapitza |Kap37| and Allen and Misener |AM37| : the discovery 
of ^He superfluidity by Osheroff et al. |()RL72[ in(;RT;72| . For later refer- 
ence, we would also like to cite the production of helium clusters by su- 
personic expansion of a helium gas by Becker et al. |BKL61| -which finally 
led to the observation of the extremelly-weakly-bound helium dimer by Luo 
et al. iLMK+QSj . and by Schollkopf and Toennies |ST94] - and the production 
of effectively low-dimensional helium systems by adsorbtion on a graphite by 
Bretz et al. |BDH+73| . 

A few theoretical works worth being cited as well are: the study of the 
role of the excitation spectrum in determining the superfluid behavior of '^He 
by Landau |La.n41[ ILan47| and the introduction of correlations to account 
for the strong repulsion at short distances between two helium atoms by 
Jastrow |.Ia.s55| and by Feynman |Fey54t lF('56| . Finally, we would also like 
to emphasize the original application of density functional techniques to the 
study of atomic quantum liquids by Stringari |Str841 lDS85| . 
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Recently, a great deal of work has been devoted to study quantum liquids 
in restricted geometries |KNfl2| . One important feature of these systems 
is that their internal structure becomes more easily observable than in bulk 
liquids due to the restricted motion of the particles in the confining potential. 
Among these systems the study of quantum films has received particular 
attention. They consist of liquid helium adsorbed to a more-or-less attractive 
fiat surface. In 1973, M. Bretz et al. |BDH^22| observed for the first time 
the adsorption of ^He onto the basal plane of graphite. In the last few years, 
adsorption properties of helium on different substrates such as carbon, alkali 
and alkaline-earth fiat surfaces, carbon nanotubes and aerogels have become 
a fertile topic of research. 

The structure and growth of thin films of ^He adsorbed to a substrate was 
studied by Clements et al. [CEKSQB] employing the optimized hypernetted- 
chain Euler-Lagrange theory with realistic atom-atom interactions. It turns 
out that films with low surface coverages (where all atoms cover the surface 
with a thickness corresponding to a single atom), can be reasonably well de- 
scribed by a two-dimensional model. In connection with these systems, an 
interesting question naturally arises as how physics depends on the dimen- 
sionality of the space. 

The homogeneous 2D liquid has been studied using different theoretical 
methods, such as molecular dynamics fCSZlJ and quantum Monte Carlo sim- 
ulations, either Green's Function ^CK8Sj or diffusion |2.BC96aJ techniques. 
The inhomogeneous case was studied by Krishnamachari and Chester who 
used a shadow variational wave function to describe 2D puddles of liquid 
^He [KC99| . 

In the following two chapters we present our studies on ^He systems in two 
dimensions (2D). In Chapter IHl we report results of Variational (VMC) and 
Diffusion Monte Carlo (dmc) calculations. First, we introduce the Monte 
Carlo techniques. Then, we show our results for 2D clusters with a finite 
number of helium atoms. We study both the energetics and the density 
profiles of these systems. 

In Chapter [3 we present a similar study developed in the framework of 
the Density Functional techniques, which allow us to study much larger clus- 
ters of ^He. We start the chapter giving a short introduction to theoretical 
background of our work: the Hohenberg-Kohn theorem and Density Func- 
tional theory. Then we present how to use the data obtained from the DMC 
calculations to set up a zero-range Density Functional. With this functional, 
we perform calculations for large clusters, which would be computationally 
prohibitive for a DMC calculation. Finally, we compare the results obtained 
with both methods. 



Chapter 6 



Quantum Monte Carlo study of 
two-dimensional ^He clusters 

Your bait of falsehood takes this carp of truth: 
And thus do we of wisdom and of reach, 
With windlasses and with assays of bias, 
By indirections find directions out. 

Williann Shakespeare, Hamlet (II, 1) 

6.1 Short introduction to Quantum Monte Carlo 
methods 

Quantum Monte Carlo (qmc) methods are powerful numerical techniques to 
solve the Schrodinger equation for interacting many-body systems. There 
is a large variety of these methods: from the easy-to-program Variational 
Monte Carlo (VMC), to the more powerful Diffusion Monte Carlo (dmc) 
and Green's Function Monte Carlo (gfmc) and the more sophisticated Path 
Integral Monte Carlo (PIMC) that is able to cope with finite-temperature 
problems. One can find in the literature a number of introductions to and 
comparisons among the various methods, see e. g. |Cep95[ IGua,98| IHLR .94J 
and also |Astfl4| rWil96| . Here we will just present a concise introduction to 
the philosophy and algorithms actually used in our calculations (vMC and 
dmc). 
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6.1.1 Importance sampling and Metropolis algorithm 

Let us denote by $ a wave function corresponding to some state of our 
interacting system. The postulates of Quantum Mechanics say that 

\^(R 

with R = {ri,r2,--- ,TAr}, is the probability that the N particles of the 
system (in our case, ^He atoms in a droplet) are located within a volum d-R 
around the positions {vi, r2, ■ ■ ■ , Vjy}, at time t |GP89j IMes99| . Indeed, $ 
contains all the information that could in principle be retrieved from the 
system. 

However, we are usually interested only on averages of operators weighed 
with p(-R), 



(A) = = J dRpiR)AiR) . 



For example, the energy of a state is just the expectation value of the Hamil- 
tonian in such a state: 

Em = ^ ' = T + V = 
L J ($1$) 

_ ElJ-dR<l>*{R)^HR) E.<JdRVir.,rmR)\' 

where we assumed that only two-body interactions V are present. In order 
to know E[^] not all the information contained in $ is required, but just 
some information on its spatial variations (to evaluate the kinetic energy 
T) and the two-point distribution function obtained once the coordinates of 
N — 2 particles have been integrated out. Therefore, two questions must be 
addressed: (1) how to find a good approximation to the wave function $ of 
interest (usually, the ground-state wave function); and (2) how to calculate 
these averages once the wave function is given. 

One could naively think in a R^'^-generalization of the simple trapezoidal 
approximation with M points to a one-dimensional integral, 



f dx/(x) ^ ^^^/(xfc) , Xk = a + ^{b- a) 
k=i 



the main difference being that now the integration is over Nd dimensions. 
However, the exponential increase in size of the (hyper)volume of integration 
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with the number of particles in the system, together with the fact that most of 
this (hyper)space usually contributes poorly to the total value of the integral, 
make such 'traditional' integration rules (trapezoidal, Simpson's, orthogonal 
polynomials, etc.) unapplicable, and one has to resort to more sophisticated 
methods such as Monte Carlo techniques. 

Here one evaluates averages such as those in (j6.2|] stochastically: a number 
Nw of points (called walkers) {Rk} {k = 1, - ■ ■ , Nw) are randomly drawn in 
R^"^ according to the probability distribution p{R). Therefore, more points 
are generated where they have more weight, and very few (if any) where the 
contribution to the integral is small. The fact that the points are sorted 
according to is what makes these quadrature techniques so powerful to 
perform multi-dimensional integrations. One usually refers to the function 
used to weigh the configurations (in the present case, $) as the 'importance 
function', and the technique 'importance sampling'. 

Then, the values of the energy of the system (or whatever the observable 
to estimate is) in such configurations = e{Rk) are recorded. For a large 
enough number of points Rk, the Central Limit theorem guarantees that 

= / dRe{R)piR) = -Lj^efc. (6.3) 

One way to obtain a set of points in M.^'^ distributed according to p{R) 
is by means of the Metropolis algorithm |MR,R,"'"53] : 

1. Start by generating N positions {vi} {i = 1, ■ ■ ■ , N) distributed ran- 
domly within configuration space. These points define an starting 
walker, R. 

2. For each position, perform a random move as 

ri = r, + , 

where £ is the typical size of the random step, and is a three- 
component vector whose elements are random numbers uniformly dis- 



its original position to a new one inside a box of side £ around it. 



tributed in [—2, g]- Thus, each particle in the walker has moved from 



The probabilities of the system being in the configuration corresponding 
to R = {ri, ■ ■ ■ , rjv} and its proposed successor R = {r[, ■ ■ ■ , r^} are 
calculated. Then, the proposed movement R ^ R' is accepted with a 
probability: 



P(R^ R') 



1 p{R') > p{R) 

p{R')/p{R) p{R')<p{R) 
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If the movement is accepted, we set = Vz. Otherwise, we keep -R 
as is; that is, the new walker is identical to its predecessor. 

4. After a certain number of iterations of steps 2-3 ('equilibration loop'), 
the set of points in the last walker should be distributed according to 
p{R). Then, it is time to start saving data to perform the average 
of Eq. (16 .^l) . To this end, one performs again a number of iterations 
of steps 2-3, but now the quantity to be averaged is stored with the 
value corresponding to R after the acceptance/rejection step ('averag- 
ing loop'). That is, if the new walker is identical to its predecessor, one 
must record the same value e(-R) saved in the previous loop. 

A chart illustrating this process is in Fig. 16.11 The size of the step £ must 
be chosen carefully: too short a step will result in a large acceptance of 
moves, but also in very correlated configurations and a long time to reach 
equilibrium. On the other hand, too long a step will probably result in a large 
number of rejected moves, and also a long time for equilibration. A good 
choice is the one that results in an acceptance around 50% — 70% |Gua,98] . 
In our calculations, £ ^ 6 A. 

It is also important to note that both the equilibration and averaging 
loops have to be repeated a number of times, called blocks. Indeed, each 
configuration produced in a step of the Metropolis algorithm is strongly cor- 
related with the previous one. Thus, a direct average as suggested in (|6.3|) 
would produce typically too small a variance. In order to have uncorrelated 
data Cfc to average, one divides the simulation in a set of blocks, each of 
which composed of a certain number of simulation steps. For each block, an 
average of the energies accumulated in it is performed, resulting in a 'block 
estimate' for the the energy Cbiock- A careful analysis on the length of the 
blocks has to be performed in order to ensure that the configurations in the 
different blocks are uncorrelated, so that the Cbiock values are effectively in- 
dependent (see |FPR91 Nig98| ). Once this has been done, one can finally give 



an estimate for the energy and the associated uncertainty, 

-j^ -'Vblocks 

Cblock 



A^blocks , , , -, 
block=l 



Ae 



\ 



ocks 

— — E ( 

A^blocks — 1 - - 



-\2 

Cbiock — ei 



block=l 



where A'biocks is the number of blocks in the averaging loop. In our simu- 
lations, typically A^biocks ~ 100, and the number of Monte Carlo steps per 

block AT^teps ~ 1000. 
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Figure 6.1: Flow chart of the Metropolis algorithm: equilibration process 
(left) and averaging process (right). 
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6.1.2 Variational Monte Carlo 

The Variational Monte Carlo (VMc) method relies on the 'Variational The- 
orem', which can be stated in this way: 

Variational theorem. The expectation value of the Hamiltonian over an 
arbitrary wave function $ is larger than, or equal to, its expectation value on 
the exact ground state, Eq. The equality only holds if ^ coincides with the 
exact ground state wave function, \l/o. 

Proof. This theorem is easily proved expanding $ on the basis of eigenstates 
of the Hamiltonian {^'n} = 0, 1, 2, ■ ■ ■ ) (with corresponding eigenenergies 
{En})- Without loss of generality, we assume that $ is already normalized 
to unity, so that we can write 

n 

c„ |c„|e[0,l], J2\cn\^^l. 

n 

Now, we insert this expression into the expectation value of the Hamiltonian: 

nm nm 

^ ^ (^nPmEm^nm — ^ ] \ ^n\ ^ -^0 ) 

nm n 

where we assumed that the eigenstates are properly orthonormalized and 
ordered in such a way that Eq < E^ < • ■ ■ . □ 

In a VMC calculation, one proposes a trial wave function $t with one 
or more (say p) variational parameters {ttj} {i = 1, • • • ,p). The functional 
form assumed should be based on all the available physical information of 
the system, and not too many variational parameters should be included in 
order to simplify the minimization process. The scope is to find the values 
of the parameters {ai} that minimize the expectation value 



Eriia,}) := E[^t] = M' , ' ^ ■ (6-4) 



($t| 


\H\^t) 


($^1 





The corresponding wave function ^({aj}) will be the best approximation to 
^0 with the functional form of $. 
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One evaluates the expectation value stochastically as outlined 

above. For the case of two-body interactions, and taking for simplicity a 
central interparticle potential, V{ri,rj) = V{\ri — r^l), we have 

^'«)-i:|;ff^^Enk.-..i). (6.5) 

1=1 ^ \ J 

Note that this expression is only meaningful for $r(-R) 7^ 0. This is not a 
problem since configurations for which = are in fact discarded by 

the Metropolis algorithm. There will be no problems in implementing this 
algorithm for Bose sytems, whose ground states are known to have no ze- 
roes |MPMfl2] . The situation is more complex for fermionic systems, where 
Pauli's principle forces "^{R) to vanish whenever R contains two indistin- 
guishable particles in the same position, e. g. $(ri, ■ ■ ■ , r, ■ ■ ■ , r, ■ ■ ■ , rjy) = 
0. This is known as the 'sign problem' for fermionic systems, and a number 
of methods have been developed to overcome it. However, as we are inter- 
ested in the study of "^He, which is bosonic and will not suffer this problem, 
we will not discuss further this point, but refer the reader to the litera- 
ture |Cep981 IHT:R94[ IWilQBL ICTm] . 



A good feature of VMC calculations is the 'control' on the wave function 
that one has. Indeed, the calculated value ()6.4^1 will depend on the functional 
form of $T and the values of its variational parameters a,. Therefore, one 
can get insight into the physics of the problem by checking different forms 
for $T, which parameters affect more the value of and, finally, what are 
the optimal values for these parameters. 

Also, if one has been able to find the exact ground-state wave function, 
\l/o5 one can make averages as that in Eq. (I6.2jl for any operator in order to 
retrieve more information on the system, which might be too difficult to ex- 
tract analytically. For example, the exact ground state of a one-dimensional 
system of hard-core bosons is known to be the absolute value of a Slater de- 
terminant of plane waves. In this case, despite having the exact "^o, analytical 
calculations are quite difficult, and Monte Carlo methods have been useful 
to obtain its momentum distribution, correlation functions, etc. [Ast04j. 

As a drawback, VMC is unable to correct by itself any 'bad' behavior 
(symmetries, long- or short- range orders, etc.) that we may have introduced 
in $. That is, if we have assumed a wave function with a different symmetry 
than ^0, the VMC calculation will be equally valid, but the upper bound 
we obtain for the ground state energy may be far from its true value. Any 
other expectation value that we calculate will be equally biased. A way to 
overcome these problems is to perform a Diffusion Monte Carlo, calculation 
which is able to provide the exact ground state energy for a Bose system. 
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6.1.3 Diffusion Monte Carlo 

The Diffusion Monte Carlo (dmc) metfiod consists of a stocliastic solution 
of the Schrodinger equation, 

d 

ih—^{R,t) = H<^{R,t) . 

We introduce an 'imaginary time' variable* r = it, so that 

t) = {H~ EtMR, t) , (6.6) 

being Et an (arbitrary) energy shift whose utility will become clear below. 
The formal time-dependent solution to this equation can be easily expressed 
in terms of the eigenfunctions and eigenenergies of the Hamiltonian as 

oo 

<^{R,t) = Y,Cn^n{R,r) = J]c„e-(^-^-)^/'^vl/„(i2,0) 

n=0 n 

= J2 c„e-(^"-^-)^/'^^„(i?, 0) ™ Co^oiR, 0)e"(^°"^-)^/^ . 

n 

Therefore, a trial function propagated on imaginary time, will converge to 
the ground state of the system, as long as it is not orthogonal to it (cq ^ 0). 
In general, the propagation will converge to the less energetic state with a 
non- vanishing overlap with while the other components are exponentially 
suppressed.^ Moreover, if the energy shift Ex equals the ground-state energy 
Eq, $(i?, r) has a stationary behavior at large r, which coincides with the 
ground state solution (up to a normalization constant). 

Let us see now how to perform this propagation on imaginary time. To 
this end, we consider a Hamiltonian with two-particle interactions V and a 
possible external potential \4xt, 

^ = E + E ^-(r.) + J2 nr.,) 

i=l i=l i<j 

where := \ri — rj\. The Schrodinger equation in imaginary time for this 
Hamiltonian can be cast in the form 

-h-^^{R,t) = -DV]i<!> + V{R)^ -Et^ , (6.7) 

*The resemblance between this variable and the time variable of finite-temperature 
Green's functions in Sect. 11.21 and Annendix 1X1 should not lead to confusion. The latter 
can be related to the temperature of the system, while t has no physical meaning in the 
Diffusion Monte Carlo formalism. 

^This fact allows for a QMC search also of excited states, see e. g. |("B88[ [rTB(^96bj . 
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where we defined D = h^/{2m) and V{R) = Kxtl^^i) + Z]i<i ^('^uOi 

is a shorthand notation for V^. Next, we define a function / as the 
product of the exact ground-state wave function, ^ and a trial one 

f{R,r) :=$T(H)^(i2,r). 

One can optimize the convergence of the DMC procedure by starting from a 
'reasonable' trial function, typically obtained in a previous VMC calculation. 
Substitution of this expression into (|6.7p gives an equation for / that is fully 
equivalent to the Schrodinger equation, 

-h^f = -DVy + DVr ■ [Ff] + (Eioe - ET)f , (6.8) 
with the definitions 

2 

F := -T—^R^T quantum drift force, (6.9a) 

Eloc := -^H^T local energy. (6.9b) 

The use of / instead of \l/ is recommended because it reduces the variance of 
the calculated observables. 

The formal solution of Eq. ()6.8^1 in i?-space can be written as 

f{R,T) = {R\f{T)) = I dR'{R\e~^^-^-^^/'^\R){R\fm ^ 

= j dR'GiR,R',T)f{R',0) , (6.10) 

where we have introduced the Green's function 

G{R,R',t) := (i?|e-(^-^^)"/^|i?') . 

Up to this point, no approximations have been made, which means that if 
we were able to know the exact G, we could get f{R, r) for any r and obtain 
from it the ground state wave function However, usually one does not 
know G exactly, but it can only be found for small time steps. Then, the 
value of /(-R, r) for r — >^ oo is found iteratively, 

f{R, t + 5t) = j dR'GiR, R', 5T)f{R', r) . 



■'■We drop the subscript '0' for notational simplicity. 
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To obtain a reasonable approximation for G, let us split the operator 
acting of / in Eq. (I6.8|l into three parts as follows 

Hi = -DV]^ (6.11a) 
H2 = D[{Vr-F) + F-Vr] (6.11b) 
Hs = Ei,,{R) - Et . (6.11c) 

We define also the Green's functions characteristic of these operators: 

GiiR,R',T) := {R\e~^^^/^\R') . 

Even though the operators Hi do not commute with each other, we can 
use the Baker-Campbell-Hausdorff formula [Clam 981 IRak02l IHau06l IWi167l 
Ryd96| , 

expXexpF = exp(^X + F + |[X,f] + ---) , (6.12) 

to approximate the exponential to some fixed order for short time lapses 6t; 
for example, to first order, we have 



(6.13) 



Similar expressions can be built with a higher degree of accuracy («. e., valid 
up to a higher power of 6t). Then, the full Green's function will read 

GiR,R',T) = jj dRidR2GiiR, Ri,r)G2{Ri, R2,r)G3iR2, R',t) . 

The solutions for the separate Green's functions are easy to find. For 
example, Hi corresponds to a purely diffusive problem, and the solution is 
well-know to be that of a random walk 



Gi{R,R',t) = {AuDt 



-3N/2 



exp 



{R - R) 

ADt 



The effect of H2 is the same as that of a classical drift force F pointing 
towards the regions where the trial wave function is maximal: 

G2{R,R',t) =6{R- R{t)) , 

i?(r) such that | • 

Finally, the Green's function corresponding to H3 is 

GsiR, R, r) = expliEr - E,,,{R))r / h]6 {R - R) . 
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The effect of Gi and G2 is implemented by proposing random movements of 
walkers of the type 

Ri ^ Ri = Ri + -\ F , 

m 

where is a 3A^-component vector whose elements are random numbers 
generated according to a Normal{0, 1) distribution. The G3 term is usually 
called the 'branching term' as it is the only contribution to G that does not 
conserve the weights of the walkers as provided by ^t- walkers with lower 
Eioc adquire larger weights, while those with larger i^ioc have smaller weights. 
This is implemented in the algorithm in the following way: 

1. The factor s = exp[{ET — Eioc{R))T/h] is computed for each walker. 

2. A number ^ is randomly generated with uniform distribution in [0,1]. 

3. The integer part of nsons = s + ,^ is calculated. 

4. The walker R is replicated risons times. (In practice, we set a maximum 
value risons = 10 in our calculations.) 

Thus, walkers with lower i^^ioc 'reproduce', while those with larger i^ioc 'die 
out'. This behavior will ultimately correct any component of the assumed $t 
orthogonal to Note that if $t happens to be an eigenstate of H, Eioc{R) = 
E\oc [cf. Eq. (I6.9bjl ]: all walkers have the same weight and no 'reproduction 
effects' take place. In the chart of Fig. 16.11 this process takes place after 
updating the positions and before storing the values for the averages. 



6.1.4 Building the trial function 

From the previous discussion, we see that a crucial ingredient in both VMC 
and DMC calculations is the definition of the trial wave function $r. Indeed, 
the 'velocity' to converge from a given $t to the exact ground state wave 
function will depend on the 'distance' between these two functions, and the 
accuracy of the average values obtained with the converged / will also depend 
on it. Therefore, we must construct $r carefully. 

For a Hamiltonian with central, two-body interactions, even if they are 
strong as in the case of helium, experience has shown that a good guess has 
the Jastrow form |.Tas55l IPM02| 



$T = A/" 



■ N 



(6.14) 
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Here Af is the normalization constant and /i is a one-body function ac- 
counting for the effects of a possible external potential. In our case of 
two-dimensional droplets, we simulate the self-trapping of the systems by 
introducing single-particle factors of the forms 



/i(r) = exp 



a 



/i(r) = exp [-ar] 



harmonic; 
exponential. 



(6.15a) 
(6.15b) 



As it will be shown later, however, only the translationally invariant part of 
n /i is to be used as there is no external trapping potential that locates the 
droplet in space. 

On the other hand, /2 contains two-body correlations and should go to 
unity for large inter-particle distances. In particular, for studies on liquid 
helium the following form introduced by McMillan |McM65j has shown to be 
useful:^ 

(6.16) 

Here h and p act as variational parameters with respect to which the expec- 
tation value {^t\H\^t) will be minimized. For helium in three dimensions, 
the values 6 ~ 3 A and z/ fa 5 provide a good Ansatz for the calculations, 
as they account for the strong repulsion at short distances r < a = 2.556 
A characteristic of the helium-helium interaction, which in our calculations 
is described by the Aziz HFD-B(HE) potential |AMW87| . see Fig. O As 
these values are related basically to the two-body potential, we think they 
will also be a reasonable starting point for our two-dimensional case. We 
note that these values give the correct behavior of the correlation function 
at short distances but they do not produce the correct behavior at long dis- 
tances. However, we expect this not to be so important for finite systems 
such as the small droplets we will be dealing with |PPW86j . 



more general approach would be that of searching the optimal functions /i and /2 
by solving the Euler-Lagrange equations d{H)/Sfi = 0, see ^M02j. We will not attempt 
this here, as the functional forms for the fi that we use have already been proved to be a 
good starting point for VMC calculations with helium and, ultimately, the dmc calculation 
in Sect. 16 .31 will 'correct' any bad behavior of the trial wave function. For the same reason, 
and in order to ease the operations to be performed on it, we are not including three-body 
correlations into <I>t. 



/2(r) 



exp 



- - 

2 \r 
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Figure 6.2: Aziz HFD-B(HE) interaction potential for ^He atoms. 
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6.2 VMC ground-state energies 

To study a system of N ^He atoms in two dimensions we start from the 
following trial wave function 



$1 



exp 



i<j 



V 2 

« 2 



(6.17) 



The first term inside the brackets corresponds to the simple McMillan two- 
body correlation, while the second term is the translationally invariant part 
of a harmonic oscillator (ho) wave function with parameter a, extracted from 
the relation 

i=l i<j 



"-CM > 



to roughly confine the system to a region of radius ~ a~^. 

In our calculations, we employ the value h'^/m4 — 12.1194 K for the 
atom mass. The parameters b and have been fixed to the values 3.00 A and 
5, respectively thus restricting the variational search to the HO parameter 
a. The optimal values found as a function of the number of particles in 
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Table 6.1: Variational results for the ground-state energy per particle E /N 
of 2D *^He puddles of various cluster sizes. The optimal confining HO param- 
eter a is given in and all energies are in K. The expectation values of the 
kinetic and the potential energies are also displayed. The column labelled 
KC refers to the VMC results of Ref |KC99| . The numbers in parenthesis 
stand for the statistical uncertainty on the last figure. 



N 


a 


E/N 


T/N 


V/N 


KC 


8 


0.1565 


-0.2239(2) 


1.3003(6) 


-1.5242(5) 




16 


0.129 


-0.3510(2) 


1.7354(6) 


-2.0864(5) 


-0.380(8) 


36 


0.094 


-0.4532(4) 


2.031(3) 


-2.484(3) 


-0.471(7) 


64 


0.073 


-0.4961(7) 


2.159(2) 


-2.655(2) 


-0.528(5) 


121 


0.054 


-0.5241(6) 


2.223(2) 


-2.747(2) 


-0.570(7) 


165 


0.047 


-0.5328(3) 


2.289(1) 


-2.822(1) 


-0.602(7) 


512 


0.0266 


-0.5493(5) 


2.282(3) 


-2.831(3) 


-0.621(2) 


oo 


0.0000 


-0.6904(8) 


4.312(2) 


-5.003(1) 





the droplet are given in Table l6T] together with the expectation value of the 
Hamiltonian and the separate contributions of kinetic and potential energies. 
It can be seen that the total energy results from an important cancellation 
between these contributions, which is in fact larger than in the 3D case. Let 
us recall that in 3D bulk, the energy per particle results from adding 14 K 
of kinetic energy to ~ —21 K of potential energy. In 2D, both kinetic and 
potential contributions are very similar, hence calculations are more delicate 
owing to the statistical uncertainties. 

The last column of Table 16.11 reports the VMC results of Krishnamachari 
and Chester |K(y99| . As compared with their results, our calculations pro- 
vide smaller binding energies in spite of the fact that the interaction used 
in |KC99j is an older version of the Aziz potential, which tends to underbind 
the systems. This is due to their use of shadow wave functions, which con- 
tain more elaborate correlations not present in our simple trial wave function. 
Finally, the VMC energy for the bulk system corresponds to the saturation 
density obtained in the DMC calculation of Ref |GBC96aj . po = 

0.04344 A-2. 

We have also performed calculations using a diff'erent trial wave function, 
replacing the translationally invariant HO part by an exponential one, i. e. 

<i>t(r) = n^^p 

to check whether a larger tail in the wave function results in more binding. 



a 



(6.18) 
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Actually, we do not find significant differences for small values of A^. For 
instance, in the case = 8, using the same values for b and u as before, 
we get E/N = -0.2178(5) K, T/N = 1.266(2) K, V/N = -1.484(2) K for 
a = 0.035 A~^. When the values of b, v and a are all optimized, we obtain a 
slightly larger binding energy, EjN = -0.2267(8) K for 6 = 3.04 A, z/ = 5.0 
and a = 0.035 A~^. For greater values of A^, the harmonic Ansatz tends 
to provide more binding than the exponential. For example, with the expo- 
nential Ansatz, for A^ = 16 we get E/N = -0.1816(7) K for a = 0.023 k-\ 
and optimizing the different parameters we get E/N = —0.2514(6) K, with 
b = 3.04 A. In conclusion, the correlated HO wave function (|6.17p seems good 
enough to be used as importance function in the DMC calculations. 

6.3 DMC ground-state energies 

After finding appropriate trial wave functions $T(a) with the previous VMC 
calculation, we will now perform DMC calculations to find the exact ground 
state energy and density profiles of the ^He two-dimensional clusters. We 
have used both first and second order |BC94j propagators [see Eq. ()6.13|) ] in 
the present work and both of them provide the same extrapolated energy, 
within statistical uncertainties, using the optimal trial function of the VMC 
calculation as guiding function. 

Our simulations have been carried out with a population of typically 400 
walkers. As usual, some 'equilibration' runs are first done to establish the 
asymptotic region of the short time propagator. Then the average blocks 
are performed. This procedure is done for several values of the time step 
6t. Finally a fit of the different obtained energies E^Mci^T), either linear 
or quadratic, has been carried out to obtain the extrapolated energy, which 
is the one we report. For example, for A^ = 16 the time steps 6t =0.0001, 
0.0002, 0.0003, and 0.0004 have been used to perform the extrapolation. In 
general, the statistical uncertainty for each time step is of the order of the 
uncertainty in the extrapolated value. 

In Table HT^ we present the results of our linear and quadratic DMC cal- 
culations of the total energy per particle for puddles containing A^ atoms. 
The linear and quadratic DMC results are compatible within their error bars. 
We have reproduced the results of the binding energy per particle of homo- 
geneous 2D liquid *^He at the equilibrium density Po"° = 0.04344(2) A'^ 
obtained in Ref. |GBC96a] . where the same version of the Aziz potential 
was used. For this case, the simulations have been carried out for a sys- 
tem of 64 atoms with periodic boundary conditions, for which the errors due 
to finite size effects are smaller than the statistical uncertainties |WCK88j . 
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Table 6.2: Energy per particle (in K) for 2D ^He puddles for various cluster 
sizes obtained with the linear and quadratic DMC algorithms. The figures in 
parenthesis stand for the statistical uncertainty. 



N 


linear 


quadratic 


8 


-0.2613(4) 


-0.2612(2) 


16 


-0.4263(4) 


-0.426(1) 


36 


-0.578(2) 


-0.575(3) 


64 


-0.658(4) 


-0.652(4) 


121 


-0.710(2) 




oo 


-0.899(2) 


-0.8971(6) 



As expected, the DMC results lower the corresponding energies obtained by 
VMC either with our simple variational wave function or with a shadow wave 
function |KC99j up to ~ 25% in the case of the bulk system. It is worth 
stressing that the final DMC result for the energy does not depend on the 
starting trial wave function for a boson system, a fact that in the present 
case has been numerically checked for the Gaussian and the exponential An- 
satze, Eqs. (I6.17|) and (I6.18jl . Indeed, for boson systems the DMC method 
provides the exact ground-state energy, within statistical uncertainties, given 
a long enough computing time. 

6.4 Energy and line tension: introducing the 
mass formula 

For a saturating self-bound system, the ground-state energy per particle can 
be expanded in a series of powers in the variable A^"^/"^, where N is the 
number of constituents and d is the dimensionality of the space. This is the 
well-known mass formula |vW35[ lBB36j . which in the present case reads 

E{N)/N = Eb + eiz + Ecz'^ + ■ ■■ (6.19) 

with z = N~^^'^. The two first coefficients of this expansion correspond to the 
bulk energy Eb and the line energy Ei (the 2D equivalent of the surface energy 
for 3D systems), out of which the line tension A is defined by 2nro\ = Ei. 
Here tq is the unit radius, that is the radius of a disk whose surface is equal 
to the inverse of the equilibrium density of the 2D bulk liquid, i. e. Poittq = 1. 
Finally, Ec is the so-called curvature energy. 



6.4 Energy and line tension: introducing the mass formula 117 



Figure 6.3: Energies per particle (in K) of A^-atom puddles as a function 
of A^~^/^, obtained from our VMC (empty squares) and linear DMC (filled 
circles) calculations. Solid lines stand for parabolic fits to the data. The 
dashed line is a straight line between the = 8 and bulk DMC values. 
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Our calculated ground-state energies (Tables IHUl and 16. 2^ are plotted in 
Fig. 16.31 as a function of N'^/"^ . (Regarding the DMC data, we just plot the 
results obtained with the linear algorithm for the sake of clarity.) One can 
see that the differences between our VMC and DMC energies increase with 
the number of atoms in the puddle. This shows that the DMC calculation 
proceeds as expected, getting rid of the 'incorrect' features of our simple 
guiding function, which becomes ever less accurate as we add more particles 
(and, therefore, more correlations) to the system. This could be improved by 
including, for example, three-body correlations in the guiding function for a 
purely VMC calculation, but nevertheless it is adequate for the importance 
sampling in the present DMC calculation. 

We have fitted these energies to a parabolic mass formula like Eq. (16.1911 . 
The coefficients of the fit are given in Table 16. 3L together with the deduced 
line tension. Notice that the coefficient Sb of the DMC calculation is identical, 
within statistical uncertainties, to the bulk energy per particle of Table [221 
on the contrary, the same does not happen for the VMC values. This gives 
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Table 6.3: Coefficients (in K) of a parabolic fit of the mass formula, as 
given in Eq. (j6.19j) . The last column displays the deduced line tension (in 

K A-i). 



Method 


£b £l £c 


A 


VMC 
DMC 


-0.654(1) 1.41(1) -0.62(2) 
-0.898(2) 2.05(2) -0.71(3) 


0.083(1) 
0.121(1) 



us confidence on Eq. (|6.19jl being suitable for the fit of the DMC data. In 
fact, the of the DMC fit is very small, = 5.7 x 10^^. Regarding the 
line tension, and despite using a different version of the Aziz potential and 
a different trial function, we observe that our VMC estimate is rather close 
to the one reported in |KC99] . Akc = 0.07 K/A. However, both VMC results 
are remarkably different from the DMC line tension that we find, namely 
Admc = (0.121 ±0.001) K/A. 

In all cases, the line energy coefficient is approximately minus twice the 
volume energy term, similarly to the 3D case |PPW86| . Also the curvature 
coefficient is not small and therefore one expects curvature effects to be im- 
portant. To stress this fact, we have also plotted in the figure a straight 
dashed line between the N = 8 and bulk DMC values. In fact, a linear fit of 
the DMC energies gives coefficients Sb = —0.885 K and ei = 1.80 K, which 
are appreciably different from the previous ones. The bulk energy extrapo- 
lated from this linear fit differs from the directly calculated value, and the 
corresponding line energy is closer to the variational one. 

We note that in both VMC and DMC cases the extracted Ec is negative, 
i. e. the binding energy is a convex function of z as it also happens for the 
3D clusters |PPW86j . This is in contrast with the value of Ec reported in 
Ref. |KC99] which was positive but rather smaller in absolute value and with 
larger error bars. Actually, as argued in Ref. |PPW86| for 3D clusters, one 
would expect the curvature correction to the energy to be positive, as it 
corresponds physically to a lack of binding of the surface atoms due to the 
curvature of the surface (see Fig. I6.4|l and, so, to an increase in the total 
energy of the system. Therefore, one should take the extracted value for 
Ec with certain caution and not emphasize its physical significance. Most 
probably it is due to the small size of the clusters used to build the mass 
formula, which ultimately forces the E/N curve to approach the abscisa axis. 
In fact, we will see in Chapter [3 that Ec turns out to be positive when the 
smaller droplets are discarded and much larger ones are taken into account. 

In any case, we have checked that, with the present DMC data, the value 
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Figure 6.4: Why do we loose binding energy in a finite system? On the left 
side we see a semiinfinite medium (shaded region) in contact with a vacuum 
(white): the lack of particles on the right-hand-side (indicated by the arrows) 
of the interface (solid line) gives rise to a reduction of their binding energy 
and, therefore, also of the total binding energy of the system. On the right 
graph we see what happens when the interface is curved (from the dashed 
line to the solid one): even more binding is lost due to the further decrease 
in the number of neighbors for the atoms close to the frontier. 




and sign of Ec are stable against different possible fits, e. g., changing the 
number of data points used to build the fit, or using a cubic mass formula. 
Also, the values of the two first coefficients eb,ei are quite robust against 
all performed fits. As an illustration, if one takes out the bulk point, the 
predicted bulk energy per particle and line tension are equal to the reported 
ones within 1% and 5% respectively. Therefore, we believe that the extracted 
line tension should be reliable, as it also happens for VMC calculations of 
three-dimensional clusters |PPW86] . where the extracted surface tension is 
in agreement with the experimental one even if < in the fit. 

6.5 VMC and DMC density profiles 
6.5.1 Calculating pure estimators 

The calculation of observables given by operators that do not commute with 
the Hamiltonian poses a new problem to the DMC method. After conver- 
gence, the walkers are distributed according to the so-called mixed probabil- 
ity distribution given by the product of the exact and the trial wave func- 
tions, f{R,T oo) = \1/(-R)<I't(-R). Therefore averaging the local values 
of the operator does not give the exact expectation value unless the opera- 
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tor commutes with the Hamiltonian. The result obtained by straightforward 
averaging is the so-called mixed estimator, 

which is of first order error in the trial wave function. A first option devised 
to overcome this problem is the use of extrapolated estimators |HLR94| . 



( yl ( -R) ) extrapolated 

2{A{R)) 

mixed {Am 



variational 



where the variational estimate is calculated with $t alone. However, also ex- 
trapolated estimates are trial- function dependent and biased |HLR94llCB 95j. 

Several alternatives have been presented in the literature in order to ob- 
tain 'unbiased' or 'pure' («. e., trial-function-independent and statistically ex- 
act) values. We have used the forward or future walking technique |HLR,94 j to 
calculate unbiased density profiles. The key ingredient to correct the mixed 
estimator is to include as a weight in the sampling the quotient \i/(R)/$T(R') 
for each walker, 

{^\A\^) _ {mj^m _ 
{^m ~ (^i^i$t) 

The value of this weight for each walker is given by its asymptotic number 
of descendants |LKC74| : 

nsons(-R;r ->CX)) . (6.20) 



$T(i?) 



Various algorithms have been proposed in order to compute this quantity. In 
this work we use an algorithm by Casulleras and Boronat |(]B95| that con- 
stitutes a simple and efficient implementation of the future walking method. 

Let us assume that we are already in the asymptotic regime where the 
usual Monte Carlo procedure would start to accumulate data to compute 
(mixed) averages. We define a new variable P that accumulates the value of 
the quantity of interest, A{R), for some Monte Carlo steps (labeled here by 
z/ = 1, ■ ■ ■ , M) within this asymptotic regime 



M 



p, = j2 Am. 



Here the same value A{Ri) is transmitted to all descendants of walker 
while the death of a walker Rj at any time implies that no more values will 
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be added to the corresponding Pj at later times. Thus, each value A{Ri) will 
appear in the set of {Pi} {i = 1, ■ ■ ■ , Nf) at the end of the interval as many 
times as descendants of the walker Ri have appeared, i. e., with a weight 
oc risons(-Rj), which is the desired one [cf. Eq. ()6.2n|) ]. Therefore, averaging 
Pi at the end of this interval we have 



This scheme is easily implemented in a usual DMC calculation, as one 
can determine a value for P at the end of each averaging block setting M = 
number of steps in a block. At the end of the simulation, from the set of 
P over the different blocks we can get an estimate for {A)p^^e and its sta- 
tistical uncertainty. A study of the dependence of P on M, typically shows 
an initial transient period during which the estimator 'gets rid' of the prop- 
erties coming from the guiding function. After this period, and for a wide 
range of block lengths M, the computed average gives an unbiased estimate 
of the expectation value of the operator |CB95j. Therefore, studying this 
dependence, one can adjust M so as to obtain a pure estimate with a small 
statistical uncertainty. 

6.5.2 Numerical results 

To obtain pure DMC estimates for the density profiles of ^He puddles, we have 
performed a careful analysis of the behavior of the profiles with varying block 
lenghts M over which they were calculated. Typically, a range M = 100 — 
1000 has been explored, showing that for M = 500 the computed quantity 
had already converged. For example, we show in Fig. 16.51 the estimated 
density profiles for block lengths M = 100 (grey line) and M = 1000 (black 
line). The improvement from one to the other is noticeable: the estimation 
obtained with a larger block length has a steeper surface, which indicates a 
larger line tension. 

The converged pure DMC estimates that we obtain in this way (with a 
larger number of blocks, to reduce the statistical uncertainties that show up 
as oscillations in the profiles in Fig. 16. 5p for the density profiles of several 
puddles are plotted with symbols in Fig. 16.61 This figure also contains an hori- 
zontal line which indicates the saturation density of the homogeneous system, 
^DMc ^ 0.04344 A"2. For the puddle containing 36 atoms, the VMC profile ob- 
tained from a Gaussian Ansatz [Eq. (I6.17jl ] is also shown as a dashed line for 
comparison. As one can appreciate in the figure, the process of optimization 
implied by the DMC method changes the profile reducing its thickness, i. e.. 




(6.21) 
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Figure 6.5: Convergence of the future walking method: density profiles for 
a puddle containing TV = 36 atoms as obtained by the future walking method 
with block lenghts M = 100 (grey line) and M = 1000 (black line). 
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producing a sharper surface. Indeed, by comparing this figure with Fig. 16.51 
one can see how the pure estimator for M = 100 is still influenced by the 
guiding function, which yields a softer boundary. This residual influence of 
vanishes completely for M > 500 for all puddles. 

A first observation on the DMC profiles is that for the smaller clusters 
the central density is below po, while for the larger values of N shown in 
the figure the central density appears to be above po, indicating a so-called 
leptodermous behavior |ST87bj . One expects that, increasing even more the 
number of particles, the central density will approach po from above, as in 
the 3D case |ST87b| ICK92] . It is worth noticing the oscillating behavior 
in the interior part of the density profile for = 121. It is difficult how- 
ever to decide whether these oscillations are genuine or are simply due to 
a poor statistics in evaluating the pure estimator at short central densities. 
Unfortunately, to discard this last option would require an exceedingly long 
computing time, within the scheme of this work. 

The solid lines plotted in Fig. 16. 61 are fits to our DMC profiles provided by 
a generalized Fermi function of the form: 

(6.22) 

The overall agreement with the numerical data is very good, and better 
than for a usual Fermi function (i. e., fixing z/ = 1). An exception is the 




p{r) 



PF 



6.5 VMC and DMC density profiles 



123 



Figure 6.6: Density profiles of ^He puddles with various number of atoms, 
N = 8 (circles), 16 (squares), 36 (diamonds), 64 (triangles up) and 121 
(triangles down), obtained from our pure estimators for the linear DMC cal- 
culations. The solid horizontal line indicates the saturation density of the 
homogeneous system. The dashed line is the VMC profile for = 36 with a 
Gaussian trial function. The figure also contains the fits to the data provided 
by a generalized Fermi function, as explained in the text. 
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droplet with 121 atoms, where the oscillating behavior of the data cannot be 
reproduced by the fitting function. However, we have already pointed out 
that the physical significance of these oscillations is not clear. 

The values of the parameters that best fit our data are given in Table lOl 
together with the thickness t and the root mean square (rms) radius cal- 
culated with them. As usual, t is defined as the distance over which the 
density falls from 90% to 10% of its central value; for a profile of the form of 
Eq. (I6.22jl this can be found analytically, and it turns out to depend only on 
c and u: 

O.i-i/i' - 1 

t = cln . (6.23) 
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Table 6.4: Parameters of a Fermi-profile fit to the density profiles. All 
lengths are in A and pp is in A~^. The parameter v is adimensional. 



N 


PF 


R 


c 




t 




8 


0.03740 


9.308 


2.156 


1.739 


8.166 


7.20 


16 


0.04204 


13.38 


2.656 


2.284 


9.580 


9.18 


36 


0.04305 


19.47 


3.104 


2.400 


11.11 


12.91 


64 


0.04386 


26.68 


3.783 


3.111 


13.09 


16.68 


121 


0.04304 


40.09 


5.566 


4.714 


18.52 


23.15 



Regarding the rms radius, defined through 

Jq°° drr^pir) 



drrp{r) 



we have checked that the value calculated within the DMC code and the one 
derived from the fit agree to better than 0.5%, except for the = 121 case, 
where the difference is 1% (most probably due to the oscillating behavior of 
the DMC profile). 

For a sharp surface (t <ti (r^)^^^), one could approximate the density pro- 
file by a step function p(r) = Pstcp0('^step — r). The rms radius corresponding 
to this profile is given by (r^)^/^ = r^tep/ "\/2, so that it grows with the number 
of particles as N^^'^. More precisely. 



That is, faster than in 3D, in which case grows as N'^^^. This behavior allows 
for an alternative determination of the saturation density by performing a fit 
of the calculated rms radii to the above relation. The value of po extracted 
from the slope of a linear fit to the mean square radii reported in Table HTH 
is 0.042 A~^, in reasonable agreement (taking into account that for these 
droplets t is not much smaller than (r^)^/^) with the determination from 
the calculation for the homogeneous system. It is interesting to note that a 
similar treatment with the fitted values for the parameter R, which intuitively 
gives the 'central' point of the surface, does not give a good estimate for po 
(0.012 A~^). In fact, R seems to grow faster than N^^'^. 

For these droplets, the thickness is continuously increasing with A^. As 
the finite value of the thickness for the semiinfinite system should define 
an asymptotic value for t, we expect that for larger puddles the thickness 
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will present a maximum and smoothly approach this asymptotic value from 
above, as happens in the 3D case [ST87b^ . To check this numerically, one 
needs the profiles of much larger puddles, which are unaffordable in a DMC 
calculation. However, it is possible to find them using a Density Functional, 
as we will show in Chapter 

Finally, we also remark the asymmetric character of the density profiles 
with respect to the point at which the density falls at half its value at the 
center of the droplet. Numerically, this behavior is refiected in the value of 
z/, which grows with and also in the increasing difference between the 
quantities R and (r^)^/^. 

6.6 Summary and conclusions 

In this chapter we have considered strictly two-dimensional systems of liq- 
uid '^He, which are of course an idealization of a real quantum film. They 
are nevertheless interesting because their study can enlighten the underlying 
structure of real quasi-2D systems. Of course, in the latter case, one has 
to take also into account the interaction with the substrate, which basically 
provides a global attractive potential. In the ideal 2D case, the suppression 
of the wave function component in the third dimension, produces an incre- 
ment of the global repulsion between atoms, resulting in a smaller binding 
energy per particle, and a decrease of the equilibrium density |AKn2| . 

We have calculated the binding energies of two-dimensional ^He clusters 
by means of a diffusion Monte Carlo method, and we have seen that they 
can be well fitted by a mass formula in powers oi z = N~^^'^. The analysis 
of the mass formula provides the value A = 0.121 K/A for the line tension, 
which significantly differs from the one obtained from a similar analysis of 
VMC data and the one previously reported in the literature |K(]99| . The 
quadratic term of the mass formula cannot be neglected and results in a 
negative value of the curvature energy, similarly to what happens in the 3D 
case [PPW86, CK92J. However, the studied clusters may be too small to give 
physical significance to this result. 

The density profiles obtained with the pure estimator have been fitted 
to a generalized Fermi function, and the behavior of the rms radius and the 
thickness, as well as the asymmetric character of the profile as a function of 
N have been discussed. 

Due to computing limitations on DMC calculations, we have restricted 
our study to relatively small clusters with N < 121 atoms. However, to fully 
understand the transition from finite clusters to the bulk medium, both for 
the energetics and the structure of the profiles, it seems interesting to have 
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results for much larger N. We will explore these cases in the next chapter 
by means of a density functional theory. 



Chapter 7 

Density Functional study of 
two-dimensional ^He 
non-homogeneous systems 

Pero jo, que sabia el cant secret de I'aigua, 

les lloances del foe, de la gleva i del vent, 

soc endinsat en obscura preso, 

vaig devallar per esglaons de pedra 

al clos recinte de llises parets 

i avanco sol a I'esglai del llarg crit 

que deia per les voltes el nneu nonn. 

Salvador Espriu, dins Final del laberint 

7.1 Introduction to Density Functional theory 

In this chapter we continue the study of two-dimensional ^He systems ini- 
tiated in Chapter IHl by analyzing the energetics and structure of the semi- 
infinite medium and slabs. Also, we will take a look at clusters larger than 
those studied by means of the DMC technique in order to clarify some of the 
points that have not been fully understood previously, e. g., the sign and 
value of the curvature energy Ec in the mass formula, or the behavior of the 
surface thickness when approaching the bulk medium. 

As it is well know, large systems are difficult to tackle with Monte Carlo 
methods as these deal with all c? x coordinates of the system components 
{d = dimensionality of the system, A^ = number of constituents of the sys- 
tem). There are several alternative approximate methods, from mean-field 
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calculations, where the many-body wave function is written as a direct prod- 
uct of single-particle wave functions, with the appropriate symmetry prop- 
erties related to the quantum-statistical character of the particles, to more 
involved techniques such as correlated basis functions (cbf) or the coupled- 
cluster method (ccm). A very good reference to have an overview of the 
power and flaws of all these techniques can be found in Ref. |FFKn2] . Here 
we will make use of the Density Functional (df) theory, which unites the 
power of computationally simple calculations with an insight to the physics 
of the problem. 

7.1.1 The Hohenberg-Kohn theorem 

The basis of the Density Functional theory relies on the famous Hohenberg- 
Kohn theorem which states that the ground state energy of a many-body 
system is determined solely by its one-body density p(r) |HK64[lFNMfl3j . In 
other words, the ground-state energy can be written as a functional of the 
density, E = E[p\. 

Unfortunately, the Hohenberg-Kohn theorem does not tell how to build 
this functional; it is just an existence theorem. We need to resort to our 
understanding of the problem to build it, which can be done in different ways. 
In the context of quantum liquids, one usually follows a 'phenomenological' 
approach: based on the available physical information on the problem, one 
defines a density functional 



through an energy density e[p\. The corresponding ground state is then 
found by functional minimization, usually restricted by the normalization 
of the density. This condition introduces the chemical potential into the 
problem as a Lagrange multiplier, namely 



Several functionals have been proposed to study helium systems, ei- 
ther pure ^He or ^He or mixtures thereof, see e. g. |Str851 IST87a[ lB.TH"'"93( 
lDT:P+95t |SzyOO| . The density functional we shall use is the simplest ver- 
sion of the zero-range functional intensively used in 3D calculations |Str85] . 
We have adjusted its parameters so as to reproduce some properties of the 
ground state of the homogeneous two-dimensional system as obtained in DMC 
calculations |GBC96a] . as well as the line tension extracted from the mass 




(7.1) 
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formula in the previous chapter. This procedure is discussed in Section 17^ 
together with the results for the slabs. Section 17.31 is devoted to the study 
of finite droplets, with special emphasis for those with a large number of 
atoms. In Section 17.41 a short discussion on hydrodynamic equilibrium in a 
droplet is presented, followed in Section 17751 bv a comparison of the obtained 
results with those of other QMC and DF calculations in 3D and 2D. Finally, 
in Section EEl the main conclusions are summarized. 

7.2 Semi-infinite system and slabs 

Density functionals to investigate surface properties of superfluid ^He were 
developed during the 1970's |ES75| . At zero temperature and in the ab- 
sence of currents, the order parameter of a bosonic system can be identified 
simply as the square root of the one-body density p{r), so it is natural to 
think of the energy of the system as a functional of the helium density. 
In analogy with the formalism of zero-range Skyrme interactions in nuclear 
physics |VB72| . Stringari proposed a zero-range density functional for non- 
homogeneous three-dimensional ^He systems |Str85j . A similar form, 



was soon used to study "^He surface properties |ST87a| and clusters |ST87bj . 

For the two-dimensional systems of our interest, we will use the same 
functional form. However, the parameters {b, c, 7} characterizing the func- 
tional for homogeneous systems have to be recalculated by requiring that 
the functional reproduces some known properties of the two-dimensional '^He 
system, such as the ground-state energy per particle (cq = —0.89706 K), and 
the saturation density (po = 0.04344 A~^), which have been determined in 
the framework of the Diffusion Monte Carlo method |(;R(^96allSMPN03j (see 
also Chapter inj. For the homogeneous system, the gradient terms in (I7.3jl 
vanish in the ground state, and one gets 



The third equation used to fix the parameters {b, c, 7} is obtained by im- 
posing the known value of the speed of sound at saturation, s = 92.8 
m/s |GBC96a] . This is related to the compressibility k of the system, and 
one can write it as 




(7.3) 



eo = bpo + cpl^"' 
= b+(l + ^)cpl. 
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Finally, the parameter d is fixed by demanding that the line tension of 
the semi-infinite system equals that obtained from DMC calculations for 2D 
clusters, Admc = 0.121 K/A (see Ch. IHlof this thesis). The behavior of the 
energy per particle as a function of density for the homogeneous system pro- 
vided by this functional reproduces very well the equation of state obtained 
from DMC calculations in a wide range of densities |(;R(^96al lNMPSn3| . As 
a first step in understanding non-homogeneous systems, we will consider the 
semi-infinite system. 

In order to determine the density profile p{r), in principle one should solve 
the Euler-Lagrange equation ()7.2|1 that results from minimizing the energy 
functional ^7.M . One should also introduce the Lagrange multiplier /i, which 
ensures the normalization of the density to the desired number of particles, 
and that is to be identified with the chemical potential of the system, 



8m 



2V2p |Vpp 



P 



P' 



+ 2bp + (2 + 7)cp^+'^ - 2dV^p = p 



(7.4) 



One should also impose the boundary conditions p(x —oo) = po and 
p'(x —>■ —oo) = 0, where the prime denotes derivative with respect to x, the 
coordinate perpendicular to the surface of the system (cf. Fig. 17. 

For the particular case of a zero-range functional ()7.3p . the line tension 
can be evaluated in a closed form, i. e., without solving for the density profile: 



\{d) 



2 / dp 

^0 



^2 



5m 



— + pd] {bp + cp'+^-p) 



1/2 



(7.5) 



Figure 7.1: Sketches of the semi-infinite medium and a slab: the x axis is 
perpendicular to the surfaces, the y axis the symmetry axis. The dashed line 
indicates the projection of the y axis on the density profile: at po/2 for the 
semi-infinite medium and at pc for the slab. 



7.2 Semi-infinite system and slabs 



131 



Then, imposing X{d) = Admc, one gets an implicit equation for d which can 
be solved numerically, thus yielding the last parameter of the functional. The 
parameters for the two-dimensional zero-range functional fixed in this way 
are listed in Table O 

Table 7.1: Parameters of the two-dimensional zero-range functional (|7.3p . 
b = -26.35 KA2 7 = 3.62 I 

c = 4.88 X 10^ kA2(i+t) d = 359 kA^ 



Once the functional is defined, one can study any two-dimensional ^He 
system. The first systems we have considered are two-dimensional slabs with 
varying central density pc = p{x = 0) (see Fig. 17. Ij) . For later reference, we 
note that a systematic study of three-dimensional "^He slabs with different 
density functionals has been presented in Ref. |SzyOO| . The Euler equation 
for the slabs is the same as for the semi-infinite system [Eq. ()7.4|] ]. The 
changes in the solution of the equation originate from the different geometry 
and boundary conditions which define the slab; in particular, the slab has 
two surfaces, while the semi-infinite medium only has one (located around 
X = 0). Translational symmetry implies that the density depends only on the 
coordinate perpendicular to the slab surface, which we call x again. Then, 
Vp = p'{x) and V^p = p"{x). In this way, the Euler equation ()7.4|1 can be 
expressed in a more convenient form, in which now p depends only on x, 

+ + 2bp + (2 + 7)cpi+T - 2dp" = p (7.6) 

The expression in parenthesis can be rewritten as a total derivative multi- 
plying it by p', 

_ ipT\ ^ ^ _ (pT \ ^ d_ /(pT \ 

VP VP PW dx \ p ) ' 

Thus, one can eliminate the second derivative of Eq. (17. 6|) by multiplying 
both sides of the equation by p' and integrating with respect to x, from the 
origin to a given value of x, 

X 

= p [p{x) - p(0)] . (7.7) 



Imposing the boundary conditions p{oo) = p'{oo) = 0, and considering that 
the slab is symmetric under inversion of the x-axis, and therefore p'(0) = 0, 




{p' 



im p 



+ hp^ + cp 



7+2 _ 



d{p 
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one finds that the chemical potential as a function of the central density of 
the slab reads 

/x = 6p, + cpr\ (7.8) 

It is worth to remind that this chemical potential is constant along the profile. 

Going back to Eq. ()7.7|) . and using again the fact that p'(0) = 0, one 
obtains 

. ^^ipgZA. (,.0) 
V 5fe + * 

This expression for p' is then used to calculate the number of atoms per unit 
length along the y axis (that is, the coverage) in a closed expression in terms 
of the density 



L Jo Jo \p'\ Jo \ hp + cp^+^ - p 

where the factor 2 takes into account the symmetry of the slab under the 
transformation x — > —x. This integral has a singularity when p ^ Pc that 
can be avoided integrating by parts. The final expression, free of numerical 
problems, reads 



N _ 4 / ^2 

rp. y/bp + cp^+^-p\l[b + c(l + 7)p1 - (£ + dp)c^il + 7)p^-i 
-4 / dp 



Also useful is the energy per unit length. 



— = 2 / dx 







_^(pO' 

8m p 



6p2 + cp2+^ + dip'f 



Using Eq. (I7.9jl and the previous definition of the coverage, one can finally 
express the energy per particle e = E/N in terms of the inverse of the 
coverage, x := L/N, 



e{x) = p{x) +Ax dp^ (^■^ + dp^ [bp + cp^+^ - p{x)] (7.10) 

Therefore, given a central density pc G [0,po], one can calculate the chemical 
potential, the coverage and the energy per particle of the corresponding slab. 
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Figure 7.2: (Top) Energy 
per particle (dashed line) 
and chemical potential (dot- 
dashed line) in kelvin for ^He 
slabs function of the in- 
verse of the coverage x (in 
A). The slope of the solid line 
is twice the line tension of the 
semi-infinite medium, which 
corresponds to the asymp- 
totic X ^ behavior of 
the energy per particle [see 
Eq. (frrni ]. (Bottom) Ra- 
tio between the central den- 
sity pc and the bulk equilib- 
rium density for "^He slabs as 
a function of the inverse of 

the coverage. "012345 

x = L/N [A] 

The energy per particle and the chemical potential for "^He slabs are re- 
ported in Fig. 17.21 as a function of the inverse of the coverage x. In the limit 
a; ^ one recovers the binding energy of the uniform system at Pq, which 
in turn coincides with the chemical potential. The energy per particle has 
a very clean linear behavior at the origin as illustrated by the solid straight 
line which provides a very good description of the energy per particle up to 
x ~ 1.5 A. The slope of this line can be analytically derived and turns out to 
be twice the line tension, which acts to reduce the binding energy per particle 
of the slab when increasing the inverse of the coverage. The derivation of this 
limiting behavior of the energy per particle can be easily obtained starting 
from Eq. (|7.1()jl , which defines the energy per particle as a function of x, by 
performing an expansion around x = 0. The leading terms result in 

e(J) = yUoo + 2Ax H (7.11) 

For values x > 1.5 A the energy per particle starts to bend horizontally 
and becomes a convex function, approaching zero very slowly. 

In contrast, the chemical potential is very flat at the origin, being deter- 
mined by the central density of the slab [cf. Eq. ()7.8|l ]. which varies only 
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slightly for small x. This can be seen in the lower panel of Fig. 17. 2L which 
displays the ratio of the central density of the slabs to the equilibrium density 
as a function of the inverse of the coverage. In agreement with the chemical 
potential, the central density is very flat at small values of x. The flatness of 
the central density and the chemical potential for small values of x, indicates 
that the slab approaches very slowly the limit of the inflnite system. 

Note that the central density of a slab can never go above the equilibrium 
density and is always a decreasing function of x. More generally, the density 
profiles of the slabs can be obtained from Eq. ()7.9|1 , which we now rewrite as 



X 



\ ax 



Pc 



-1 



Pc 



dp 



p{x) 



8m p 



+ d 



hp^ + cp^^"^ — pp 



1/2 



(7.12) 



valid for x > 0. The profile for x < is found by symmetry. Notice the 
presence again of a divergence when p ^ pc which can also be avoided 
integrating by parts. The profiles calculated for various central densities are 
plotted in Fig. O 

The size of the slab increases with the central density. A measure of this 
size is given by the radius R1/2, defined as the distance from the origin to the 



Figure 7.3: Density profiles for ^He slabs with central densities pc/ po=0.05, 
0.1, 0.25, 0.5, 0.75, 0.9, 0.95, where po is the equilibrium density of the bulk, 
indicated by the horizontal dotted line. 
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Figure 7.4: Radius (left panel) and thickness (right panel) of "^He slabs as 
a function of x. The symbols are the calculated data while the lines are cubic 
splines to guide the eye. 
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point where the density has decreased to half its central value. This quantity 
is shown on the left panel of Fig. 17.41 as a function of x. As expected, it 
diverges in the limit pc — > po and is mainly a decreasing function of x. 
However, it also presents a very shallow minimum around x ~ 3.3 A, after 
which it grows slightly. Physically, one can understand this last feature from 
the fact that for x > 3 A the density of the system is everywhere small (cf. 
Figs. I7.2M7.3|1 : each particle will have very few others to attach to, and the 
slab will become more extended. 

The profiles are also characterized by the thickness t, defined as the dis- 
tance between the points where the density has decreased from 90% to 10% 
of its central value. The dependence of the thickness on x is shown on the 
right part of Fig. 17.41 Up to x ~ 1, the thickness is a pretty fiat function, 
indicating that the surface of the slab is very much the same, and is actually 
the radius of the slab that grows very fast and diverges when x 0. The 
thickness is large for the slabs with larger x again due to the low density of 
these systems. When the central density is increased (x decreases), the thick- 
ness decreases until it has a minimum at x 2.75 A, with t ^ 8.6 A. Then, 
it increases again approaching a finite value at the origin, corresponding to 
the semi-infinite medium, t ~ 11.03 A. It is interesting to note the plateau 
that appears for x < 1, which may be related to the analogous behavior for 
the central density (cf. Fig. I7.2|l : the slabs are in many aspects similar to 
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the semi-infinite system already for x ^ 1. However, other properties (most 
importantly, the energy per particle) do need to go to the limit x — > to get 
rid of its finite-size dependencies. On the other hand, for x > 3 the slabs 
seem to change from a strongly correlated behavior over to a more 'dilute' 
one. 

7.3 Drops and line tension 

As a next step, we consider drops of a fixed number of atoms N. These 
where already studied by DMC techniques in Chapter El However, computa- 
tional limitations allowed us to study only small values of N. Here, we will 
take advantage of the computational feasibility of the Density Functional 
calculations and will extend the analysis to much greater N. In this way it 
will be possible to study the asymptotic behavior of several quantities which 
characterize the drops. 

Equation ()7.4|) for the profile can be rewriten in the form of a Schrodinger- 
like equation for p: 

+ 26p2 + (2 + 7)cp2+T - 2dpVV = fJ'P (7.13) 

As the number of atoms in a droplet is a well defined iV, one would need 
to be very careful in determining the central density p(0) so that the chemical 
potential adjusted exactly to N. However, it turns out that this equation is 
more efficiently solved by means of the steepest descent method jPFKWSTI] . 
An initial trial Pinput(^) is projected onto the mininum of the functional by 
propagating it in imaginary time r. In practice, one chooses a small time 
step At and iterates the equation 

p(r, r + At) ^ p(r, r) - At np{r, r) (7.14) 
p(r, 0) = Pinput(r) , 

normalizing p to the total number of atoms at each iteration. The time 
step At that governs the rate of convergence should be taken appropriately 
small in such a way that Eq. ()7.14p is a valid approximation. Convergence is 
reached when the local chemical potential — defined as the value of [Hp]/p — 
has a constant value independent of position. In our calculations, typical 
values have been At ~ 10~^, and some 10^ iterations have been necessary to 
reach convergence. This large number of iterations does not translate into 
long runs, as the local nature of the functional allows for fast calculations. 
Nevertheless, it is also possible that including in Eq. ()7.14p second or higher 
orders in At, the computation might speed up, which we have not attempted. 
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Figure 7.5: Energy per particle (empty circles) and chemical potential (full 
circles) of '^He droplets as a function of z = A^~^/^. Also shown is a quadratic 
fit (see text) of the results with > 512 (solid line). The straight short- 
dashed line is obtained when the mass formula (with terms up to z^) is used 
to calculate the chemical potential. The empty squares are the DMC energies 
from Chapter El and the dotted line corresponds to the quadratic fit to these 
results reported there. 




The energy per particle (empty circles) and the chemical potential (full 
circles) calculated for drops with N = 16, 36, 64, 121, 512, 1024, 2 500 and 
10 000 atoms are reported in Fig. 17.51 a.s a function of A^~^/^. Also shown 
for comparison are the DMC energies (empty squares) and their quadratic fit 
to a mass formula calculated in Chapter El The calculated energies of the 
droplets can again be represented very accurately with a mass formula of the 
type [cf. Eq. §M] 

e{N) =eb + eiz + ecz^ + ■■■ {z = N-^l^). (7.15) 

Contrary to the DMC calculations, where the largest droplet that we stud- 
ied had 121 atoms, here we have considered puddles with up to 10 000 atoms. 
In this way we can accurately study the behavior of the energy per particle 
for small values of z. Doing a quadratic fit to the calculated energies per 
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particle for > 512 and including also the bulk binding energy for z = 
one gets 

y = -0.897 + 2.0587Z + 0.664662^ (7.16) 

which is plotted by a solid line in the figure. One sees that Eb accurately 
reproduces the bulk energy per particle, which was used to fix the parameters 
of the functional; we note that a similar agreement was found for the DMC 
fit parameters (see Table l673|l . The value of = 2.0587 (cf. ei = 2.05 K for 
the DMC fit) corresponds to a line tension A = 0.121 K/A, which is the same 
as the value of the line tension of the semi-infinite system, used to build the 
density functional. Note, however, that the fit, even if it has been calculated 
with the data for > 512, is rather accurate down to = 36. Obviously, 
one can not expect a good agreement for = 16, which is a system where 
finite-size effects are very important. 

The linear behavior of the chemical potential as a function of z contrasts 
with the plateau that was found for the slabs. This linearity is easy to 
understand using the mass formula (|7.15p and the thermodynamic definition 
of the chemical potential /i = dE/dN, where E is the total energy of the 
droplet. With this prescription, the slope of the chemical potential as a 
function of z results to be 5//2. This behavior is illustrated by the short- 
dashed line in the figure, which nicely follows the calculated data down to 
A^ = 36. A similar analysis for the three-dimensional case, would provide a 
behavior of the chemical potential as a function of z = A^~^/^ dominated also 
by a linear component, but with a slope at the origin given by 25^/3, where 
Eg is the surface energy associated to three-dimensional clusters |ST87b| . 

It is interesting to remark that the coefficient of is positive. This 
sign corresponds to the expected loose of binding energy associated to the 
curvature of the contour of the cluster (cf. Fig. l6.4|) . and it is in contrast with 
the value of Ec obtained by fitting the DMC results (see Table I6.3|l . However, 
in that case the number of particles in the clusters used to build the fit was 
much smaller, being A^ = 121 the largest number of particles and going down 
to N=8 for the smallest one. In the present fit we have explicitly avoided the 
clusters with a small number of particles which can easily distort the results 
and we have considered only the cases with A^ > 512. 

Next thing to analyze are the density profiles, which we report in Fig. 17.61 
for different numbers of atoms. Contrary to the slabs, the central density of 
the droplets can be higher than the saturation density, which is indicated in 
the figure by an horizontal line. The profiles are well adjusted by a generalized 
Fermi function. 
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Figure 7.6: Density profiles for ^He droplets for N=16, 64, 121, 512, 1 024, 
2 500 and 10 000 atoms as fitted to a generalized Fermi profiles [see Eq. ()7.17|1 
and Table 1712] . The dotted horizontal line indicates the equilibrium density 
of the homogeneous system Pq. 
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which has an associated central density p(0) = Pf/[1 + e~^^^Y . The param- 
eters fitting the profiles for the different clusters are quoted in Table I7.2[ 
together with the thickness and root-mean-square radius obtained thereof. 

The left panel of Fig. 17.71 reports the central density of the different 
droplets as a function of z. For large values of N, p(0) is larger than the 
saturation density, i. e. the central part of the droplet is more compressed 
than the bulk system, which is sometimes referred to as a leptodermous be- 
havior |ST87b] . Of course, for — oo the central density tends to po- 
First, p(0) grows almost linearly with z, reaches a maximum, for z ~ 0.13 
(A^ ^ 60), which would correspond to the most compressed droplet and, 
finally, for N < 25, the central region of the droplets becomes rapidly less 
compressed than the bulk system. 

The root-mean-square (rms) radius corresponding to the profiles is shown 
in the right panel of Fig. 17.71 as a function of N^^'^. The expected linear 
behavior, associated to a constant average density [cf. Eq. (16.24 



2\l/2 



(7.18) 
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Table 7.2: Parameters of a generalized Fermi-profile fit [Eq. (j7.17|) ] to the 
density profiles obtained with the zero-range density functional. All lengths 
are in A and pi? is in A~^. The parameter v is adimensional. 
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36 


0.04494 


19.1852 


3.22054 


2.37516 
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64 


0.04974 


24.6826 


3.16845 


2.37072 


11.364 


16.276 


121 


0.04441 


32.8314 


3.12302 


2.33636 


11.226 


21.721 


512 


0.04392 


64.4245 


3.08867 


2.32368 


11.112 


43.533 


1024 


0.04378 


89.8497 


3.08584 


2.32997 


11.097 


61.342 


2 500 


0.04366 


138.628 


3.08552 


2.33826 


11.090 


95.678 


10 000 


0.04355 


274.026 


3.08653 


2.34670 


11.088 


191.276 



Figure 7.7: (Left) Central density of ^He droplets as a function of z = 
jY-i/2 rpj^g empty circles correspond to the results obtained with the zero- 
range density functional, while the line stands for a cubic spline fit to these 
data. The dotted horizontal line indicates the saturation density po- (Right) 
Root-mean-square radius of ^He droplets as a function of N^l'^ . The solid 
line is a linear fit to the data without independent term. The empty circles 
correspond to the results of the zero-range density functional. 
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is rather apparent. A linear fit to the calculated values (without independent 
term), from = 16 to = 10 000, provides Pq = 0.0434 A~^, in very good 
agreement with the bulk equilibrium density used to fix the parameters of 
the density functional. 



7.4 Hydro dynamic equilibrium in a droplet 

We proceed now to analyze the hydrodynamic equilibrium in these helium 
drops. To this end, it may be interesting to analyze the pressure in one 
of these clusters. The study of fluid interfaces goes back to Gibbs | Gib6l |. 
However, a clear definition of local thermodynamic variables has not yet 
been reached. It has been established that, across an interface, some vari- 
ables (such as the chemical potential) must remain constant, while for oth- 
ers (such as the pressure) this is not so, and an unambiguous definition is 
lacking for the latter, even though there have been several attempts in re- 
cent years |TGW+84[ IrI^^ mow94'l [77R971 IMRLQTI lORvM+99l [hT02] . The 



usual thermodynamic expression for the pressure at zero temperature is 

fdE\ 2 fde 



where we also used the definition of the chemical potential, /i = {dE/dN)v, 
e stands for the energy per particle and /x is the chemical potential, constant 
throughout the system. This expression depends only on intensive variables, 
and is therefore appropriate to introduce the concept of local pressure in a 
finite system. 

For our simple zero-range density functional ()7.3jl , the energy per particle 
can be written as a function of the density thus 



which, combined with the equation above, will define a 'pressure profile' p(r). 
Moreover, all the contributions to e(p) as well as pp are quantities computed 
during the imaginary-time evolution towards the ground state. Therefore, 
the evaluation of p(r) does not require any further computational effort. 

The calculated p{r) for two sample droplets (A^ = 36, 10 000) is shown 
as a function of a rescaled radial distance r/_R in Fig. 17. 8| with R the fitting 
parameter of Table 17.21 The value p = outside the drop is due to the fact 
that we are working in a vacuum at zero temperature, and therefore there 
is no vapor pressure exerted on the drop. The positive value of p in the 
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Figure 7.8: Pressure profile in the clusters with = 36 (left) and 
= 10 000 (right) atoms: the solid lines are the density profiles (left or- 
dinate axes) while the dashed lines are the calculated pressure profiles (right 
ordinate axis). The thin horizontal lines indicate the saturation density po 
and the dotted lines mean zero pressure. 
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innermost region of the droplet is related to the density there being larger 
than the saturation density. This behavior is also expected from the Young- 
Laplace law for a liquid-vapor interface, that relates the pressures in the bulk 
of each phase with the surface tension A and the mean curvature radius -Rcurv 
of the interface [MBL97^ .Saf03| . For the case of cylindrical symmetry, 

Pliq = Pvap + • (7.19) 

As mentioned above, in the zero temperature case Pvap = 0. The fact that the 
regions where p > and p > po do not coincide exactly is due to the gradient 
terms of the Density Functional, that are not present for the homogeneous 
case that fixes po. Using now Eq. ()7.19|) with A = Admc and the calculated 
values p(0) = 7.82 mK/A^ (A^ = 36) and p(0) = 0.447 mK/A^ (A^ = 10 000) 
as the values of the pressure in the liquid phase, we find the values for the 
curvature radii of the clusters with A^ = 36 and 10 000 atoms to be -Rcurv = 
15.5 A and -Rcurv = 271 A, which are rather close to the -R-values from 
the Fermi-profile fit (see Table I7.2jl . This can be seen as an indication of 
the validity of our definition for p(r), specially for large drops, where the 
hypothesis of a sharp interface that resides in the Young-Laplace law is better. 
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Figure 7.9: (Left) F(r) stands for the effective force that keeps the drop 
bound. (Right) The density profile (black line) and the effective force F(r) 
(red line) calculated to keep the drop with N = 10 000 atoms in equilibrium, 
see Eq. ^J^. 




Close to the surface, we find a negative value for the pressure, that is, 
there is a tension located near the free surface of the puddle. This is es- 
pecially clear for the larger drop, where the region where p < is closely 
restricted around the surface (note the difference in scales on the x-axes of 
both figures), while for the = 36 system the separation between a bulk and 
a surface regions is not easy to ascertain. Similar results have been found in 
various Molecular Dynamics calculations for liquid- vapor interfaces of differ- 
ent geometries (planar, cylindrical), even if using a variety of definitions for 
the local pressure (see, e. g., |TGW"'"84[ [MBL97J) . 

This picture of the pressure profile can be made more physically appealing 
by interpreting it in terms of hydrodynamic equilibrium. To this end, let us 
analyze the forces exerted on an area of a drop as depicted on the left panel 
of Fig. 17.91 The area will be in equilibrium if there exists a force F that 
balances the pressure difference on the two sides. This force must satisfy 

b (r) = — r dq) Vpyr H ) — -pyr ) = — r a0— dr = —dA— , (7.20) 

2 2 dr dr 

where dA is the surface of the shaded area, and the sign convention is such 
that F > pulls atoms to the center of the puddle. This force is ultimately 
determined by the He-He interaction. 

Thus, the pressure profile can be mapped onto a central potential for the 
helium atoms. The force calculated in this way for the drop with = 10 000 
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atoms is plotted on the right panel of Fig. I7.9I It vanishes for positions 
r < 0.94/2, corresponding to the region where p ~ Pq. in this region, the 
system is locally equivalent to the bulk, and no net forces are present. The 
force also vanishes far from the cluster: one could add a new '^He atom there 
which would not feel the presence of the drop. In the region close to the 
surface, F{r) presents a maximum at rmax ~ 0.98/2 ~ 269 A, a node at 
'"node ~ 0.99/2 ~ 272 A and a minimum at rmin R k, 275 A. The position 
of the node corresponds to the minimum in the curve for p{r), and indicates 
where an atom might be added and stay at rest. On the contrary, if it were 
added at smaller or larger r it would be displaced towards rnode- The linear 
behavior of F{r) close to rnode shows that this motion would be harmonic, as 
expected around a minimum of a potential. Note that the distances between 
these three positions are barely equal to the characteristic length of the He-He 
interaction, a = 2.556 A. 

7.5 Discussion of the results 

Let us finally make a comparison between the results obtained for clusters 
with the present DF calculations and the DMC ones in the previous chapter, 
and also between our results and those for three-dimensional slabs obtained 
from various density functionals |SzyOO| . 

7.5.1 2D and 3D slabs in Density Functional theory 

First of all, we compare our DF results for two-dimensional ^He slabs (Sect. 17^ 
with those found in the literature for three-dimensional slabs. In partic- 
ular, we will look at Ref. |SzyOO| which reports a comparison of calcula- 
tions performed with a variety of density functionals. Indeed, the functional 
e[p] in Eq. (|7.1I1 is in principle rather arbitrary, and only a comparison of 
the obtained results from different choices with experimental data (or exact 
calculations as provided, e. g., by the QMC method) can decide what func- 
tional dependence is better to reproduce certain properties. Thus, for 3D 
helium systems a 'fauna' of density functionals has developed, starting from 
the simple zero-range functional originally proposed by Stringari |Str84] . to 
more elaborate finite-range functionals such as those referred to as Orsay- 
Paris tDHPT90j, Catalonia pJH+ 93j. or Orsay-Trento |DT:P+95j in |SzyOO| . 

Of course, every functional is expected to give the same results for the 
homogeneous system, and other properties used to fix its parameters. The 
differences appear when calculating features not directly incorporated in their 
construction, such as surface properties of slabs or finite clusters. However, 
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they tend to present some general common behaviors. For example, the en- 
ergy per particle and chemical potential of 3D '^He slabs depend on the inverse 
coverage in a way quite similar to what we find in 2D systems (compare our 
Fig. O with Fig. 2 in |SzyOO| ). 

It is worth noticing, however, some intrinsic diff'erences between the 2D 
and 3D homogeneous systems. To this end, let us plot the reduced energy per 
particle e(p)/eo as a function of the reduced density p/po, where Cq and po are 
the energy particle and the density, respectively, at saturation (see Fig. lT.lfljl . 
We recall that ef^ = -7.17 K and pl^ = 0.021837 A'^. Even though the 
values for E/N are quite similar, the curvature of both lines is different, as 
can be seen from their extrapolations to higher and lower densities using a 
fit of the form 

e(p) =bp + cp^+^ . (7.21) 

This feature is reflected in the speed of sound at saturation density, which in 
2D is smaller than in 3D (s^^ = 92.8 m/s ws. s^^ = 237m/s). These different 
behaviors may be expected to affect also the structure of finite systems, where 
compressibility plays an important role in determining density profiles and 
other properties. It is also remarkable the broader range of densities where 
helium remains liquid in 2D as compared to the 3D case, cf. Fig. I7.1()[ 

The thickness of slabs in 3D is quite fiat as a function of the inverse of 
the coverage x = L/N (cf. Fig. 4 in |SzyOO| ). This behavior contrasts with 
the 2D case where we observe a minimum in the t{x) function (cf. Fig. 17. 4|) . 
This discrepancy might just be due to a broader range of inverse converages 
explored in our work as compared with |SzyOO| . In fact, our results are also 
quite fiat in the region x < 2, which relates to the larger slabs we have 
studied, Pc/po ^ 0.7. Indeed, the calculations presented in |SzyOO| are for 
slabs with pd pfP ^ 0.6. However, this does not rule out an infiuence of 
dimensionality on the behavior of the thickness, as is pointed out by the 
different values for the semi-infinite system: tgjj^j ^ 11 k vs. tg^j 7 A. 

7.5.2 2D clusters according to DF and QMC calculations 

Let us now turn to finite clusters and compare the results we have obtained 
by the DMC (Ch. ^ and DF (this chapter) methods. First, let us look at 
the energy per particle, plotted in Fig. 17.51 We see that both calculations 
give very similar values for = 121, while for < 121 the DF energies are 
always less bound than the DMC ones. On the other hand, there is a nice 
agreement of the mass formulas for > 121 {z < 0.1). This behavior is to be 
expected as the parameters defining the functional (6, c, 7, d) have been fixed 
so as to reproduce properties of the homogeneous and semi-infinite systems 
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Figure 7.10: Comparison of the 2D and 3D equations of state. The thick 
lines are DMC results, while the thin lines are their extrapolations ()7.21ll to 
higher and lower densities. The triangles and circles denote the spinodal and 
freezing points, respectively. The black curves and empty symbols refer to 
the 3D data, while the green curves and filled symbols stand for the 2D data. 
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(see Sect. 17. 2( ). For the same reason, the deviations for small clusters, where 
surface effects are most important, are not surprising. Indeed, from Fig. 17. 6| 
we see that the density profile of the droplets does not present a truly 'bulk' 
region in its interior for clusters with less that 121 atoms. Nevertheless, the 
maximum relative deviation of the energy per particle is just around 15%. 

In Fig. 17.111 we show the density profiles obtained in the VMC, DMC and 
DF approaches for the = 64 drop, for which the energy calculated by the 
last two techniques is quite similar and the fluctuations in the DMC profile 
are small. In this comparison one can see that the DF and DMC profiles 
are quite close to each other, with the surface location and thickness very 
similar as compared with the VMC ones. Thus, we can conclude that the 
description given by the zero-range Density Functional (|7.3|) is preferable 
over that of the simple VMC calculation. In this sense, it is worth mentioning 
that density profiles of 3D drops calculated with DMC and with finite-range 
density functionals show a similar agreement iBGH+flBj . Indeed, the values 
for R and (r^)^/^ of the generalized Fermi-profile fits for the DMC and DF 
data are in general in good agreement (cf. Tables lOl and rr2|l . Even though 
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Figure 7.11: Density profiles of tlie N = 64 cluster from the VMC (solid 
line), DMC (dashed line) and DF (dash-dotted line) calculations. The thin 
horizontal line is the bulk saturation density. 
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a closer look at the data shows that R in the DMC calculations is typically 
smaller than that from the DF ones, the behavior of the rms radius is very 
much the same, with differences at most of 6%. 

Regarding the thickness, one observes that it is quite constant for all DF 
profiles, with a slowly decreasing behavior towards the value for the semi- 
infinite system, igemi — 11-03 A. This is in clear contrast with what happens 
with the DMC data, which are growing with A^. In fact, the DF thickness 
turns out to be systematically smaller than the DMC one, i. e., the DF profiles 
present sharper surfaces than the DMC ones. 



7.6 Summary and conclusions 

We have constructed a density functional suitable to study non-homogeneous 
two-dimensional ^He systems. The parameters of the functional have been 
fixed by demanding it to reproduce some known properties of the homoge- 
neous and semi-infinite systems. Then, we have considered two-dimensional 
slabs and analytically shown that the line tension extracted from a mass for- 
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mula adapted to this type of geometry, is consistent with the value of the 
line tension of the semi-infinite system. We also observe that, as in three- 
dimensional slabs, the thickness and the central density of the slabs approach 
from below the values corresponding to the semi-infinite system, while the 
radius diverges. 

We have also studied the energetics and structure of 2D clusters. In 
particular, clusters with a very large number of atoms have been used to 
study the behavior of the mass formula and to establish how the system 
approaches the bulk limit. The central density of the clusters when N ^ oo 
approaches the saturation density from above and, therefore, the internal 
regions of the clusters are more compressed than the bulk system while the 
external regions have densities which would correspond to negative pressures 
or even below the spinodal point for a uniform system . This behavior reflects 
the existence of a tension on the boundary of the drop that keeps it bound. 

The profiles of the clusters are very well fitted by a generalized Fermi 
function. The thickness of the cluster slowly approaches the thickness of the 
semi-infinite system, as it also happens in the case of the slab geometry, but 
in this case from above. This behavior contrasts with the DMC data, which 
were however obtained for quite small clusters, far from the bulk limit. 

Finally, we have analyzed the linear behavior of the rms radius of the 
droplets in terms of N^^'^, and recovered the saturation density from the 
slope of this fit. We have also seen that the rms radii obtained from the DF 
calculation are in very good agreement with the DMC results. 

Altogether, we conclude that the proposed density functional is useful 
to perform reliable and computationally-affordable calculations for a variety 
of two-dimensional '^He systems and, in particular, for clusters where the 
number of particles is prohibitive for a Monte Carlo calculation. 



Conclusions and future 
perspectives 
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The purpose of this thesis has been the study of various atomic systems 
where the quantum effects are most visible, by means of a number of tools of 
many-body physics. In Chapter [l] we have presented in a general form the 
BCS theory of superconductivity for fermionic systems. We have made use 
of the well-known Green's functions' formalism to write the gap and number 
equations in a general form, capable of being used at zero and at finite 
temperature, both for density-symmetric and asymmetric systems. We have 
recovered the well-known result that the zero-temperature s-wave gap for the 
symmetric system in the weak-coupling limit has an exponential dependence 
on the interaction and the Fermi momentum, namely 



The critical temperature for the superfiuid transition has a similar behavior 
on the parameters of the system. One can write the corresponding result as 




(a < 0) . 



ksTc = 




'sym 



0.567A, 



■sym ■ 



In Chapter [2l we have analyzed the prospect for s-wave pairing in the 
BCS framework for large density asymmetries. The main conclusions of this 
chapter are: 
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The maximum asymmetry for which BCS pairing is possible is 

, BCS _ P1 ~ Pi 3 Asym 



PT + Pi 4 /i 



which in the low-density limit is a very small quantity. We have checked 
numerically this prediction by solving self-consistently the gap and 
number equations; the agreement, as shown in Fig. 12. 6| is remarkable 
for all a. 

We have also checked numerically the dependence of the value of the 
gap with temperature for the symmetric system, and seen that it follows 
the expected limiting behaviors at low temperatures T <C Tc and close 
to the critical region T < T^, cf. Fig. 12.41 

Finally, we have analyzed the prospect for a BCS transition in p-wave 
between identical fermions, mediated by another fermionic species. For 
a density-symmetric system, the gap 



Ap-wave 

— oc exp 



\2kFa J 



is very small, but it may be optimized by appropriately adjusting the 
density asymmetry. In practice, a sharp maximum appears for Oopt ~ 
0.478. 

In Chapter [3l we have studied two structures for the ground state alter- 
native to the usual BCS one: the plane-wave loff and the DFS phases. 

• As it is well-known, the LOFF phase can sustain pairing for larger asym- 
metries and has lower free energy than BCS. We have quantified this 
effect for the case of ^Li at ultralow temperatures and in the weak- 
coupling regime: we have shown that the LOFF phase if preferable to 
the BCS one for asymmetries 4% < a < 5.7%. 

• The DFS phase consists in an ellipsoidal deformation of the Fermi sur- 
faces of the two pairing species. More precisely, the surface of the most 
populated species becomes prolate {i. e., cigar-like) along a symmetry- 
breaking axis randomly chosen by the system. The surface of the less 
populated species becomes oblate (i. e., pancake-like), so that the two 
Fermi surfaces have, in the end, approached each other in the plane per- 
pendicular to the symmetry-breaking axis. In this case, only rotational 
symmetry is broken. 
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• We have numerically evaluated the gap and free energy of the DFS 
phase, and we have shown that it becomes the ground state for the ^Li 
system under study for asymmetries a > 4%. 

• We have proposed an easy scheme of detecting the DFS phase by means 
of a time-of-flight determination of the momentum distribution for an 
atomic gas confined in a spherical trap. 

It is our purpose to analyze in the future the excitation spectrum of the 
DPS state, to identify the new Goldstone mode that shold appear due to 
the breaking of rotational symmetry. One should also remember that more 
general structures for the LOFF phase are possible, which could qualitatively 
modify our conclusions. 

In Chapter 131 we have considered the effects of the presence of a bosonic 
species in a single-component fermionic system. 

• For the three-dimensional case, the induced pairing is optimized for a 
bosonic density 

2/3 

PB - 

2.88(1bb 

where Qbb is the s-wave scattering length for boson-boson collisions. 

• For the two-dimensional case, the logarithmic dependence of the transi- 
tion matrix on the collision energy requires a more delicate treatment, 
and the value of the optimal densities depends on the parameters de- 
termining the collision properties among bosons. Ebb, and between 
bosons and fermions, Ebf- We have shown that, for the case that 
both species have similar masses. Ebb plays only a minor role, and the 
optimal density ratio is given by 

-0.59 [InC- 0.35] , 

mp Bp 
niB Ebb ' 

Under these conditions, energy gaps of the order of the fermion chemical 
potential may be achieved. 

• For the case of trapped, quasi-two-dimensional atomic systems, this 
regime seems to be achievable for mixtures with attractive interspecies 
interaction qbf < such as ^^Rb-*^°K. On the other hand, for gbf > 
(as in the ^Li-^Li mixtures) one should also analyze the stability of the 
system. 
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In Chapter Hlwe have studied the real-time evolution of a Bose-Einstein 
condensate with spin degree of freedom. We have recalled that, for the case 
of spin / = 1, two possible ground-state structures are possible 

• 'Ferromagnetic': when the coupling constant in the collision channel 
with total spin F = 2 is attractive, C2 < 0, the spin tends to be locally 
maximized. This translates into a mixture of all spin components in 
each point. 

• 'Antiferromagnetic': when C2 > 0, the ground state is realized by min- 
imizing the expectation value of the spin. As a consequence, atoms 
with spin projection m = ±1 will tend to repel m = atoms. 

The dynamical evolution in the first case, which is realized for ^^Rb atoms 
in their electronic ground state, presents these characteristics: 

• It shows oscillations in the populations of the various hyperfine states 
m = 1,0, —1. At zero temperature, and for an initial state with zero 
magnetizaton, these oscillations are around the ground-state configu- 
ration (25%, 50%, 25%). At a finite temperature T = 0.2Tc, which we 
simulate via the inclusion of phase fluctuations, all components are 
equally populated after a fast damping process. 

• The density profiles feature small spin domains of a characteristic size 
^dom ~ 10 /im, that appear after ~100 ms. This domain formation may 
be interpreted as the consequence of a dynamical instability of the spin 
excitation modes of the 'ferromagnetic' system. At zero temperature, 
the domains formed by m = ±1 atoms overlap, thus conserving locally 
the initial magnetization. At finite temperature, the symmetry between 
m = 1 and m = — 1 is broken by the fluctuating phases, and different, 
interpenetrating domains of all components appear. In this case, the 
formation of domains is much faster: it occurs for times ~10 ms. 

We have also analyzed with 'antiferromagnetic' interactions for which: 

• The population oscillations are almost coherent at zero temperature, 
but present a very small relaxing behavior at finite temperature. 

• The density profiles barely evolve at T = 0.2Tc, while at zero temper- 
ature Josephson-like oscillations may be identified. 

For the future, we plan studies on the role of the initial conditions and 
on effects of temperature and external magnetic fields. Also the study of 
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condensates with higher spin values (in particular, / = 2) is interesting in 
order to interpret a number of experimental results. 

In the second part of the thesis, we have presented two microscopic stud- 
ies on ^He systems in two dimensions. In Chapter [6l we have performed 
both Variational and Diffusion Monte Carlo calculations for finite clusters. 
First we have built a good variational function by including Jastrow-like cor- 
relations together with a Gaussian term to take into account the self-binding 
of the drops. This function has then been used as importance function for 
the DMC calculation. The main conclusions that we get from these studies 
are: 

• From the results for the energy per particle as a function of the number 
of atoms in the drop, we have constructed a mass formula, 

E El Ec 

The first coefficient coincides with the calculated value for the energy 
per particle at saturation of the homogeneous system. From Ei we have 
extracted a value for the line energy due to the finite character of these 
systems, 

Admc = 0.121 K/A, 
which is the principal result of this chapter. 

• We have also analyzed the density profiles obtained by the forward 
walking method to evaluate pure estimators. The profiles can be very 
well fitted to a generalized Fermi function, and provide an alternative 
determination for the bulk saturation density, which roughly coincides 
with the value from DMC calculations. 

In Chapter [71 we have used the previous results for the homogenous 
system and the line tension to build a zero-range density functional, with an 
energy density of the form 

The parameters {b, c, 7} have been fixed so as to reproduce the energy per 
particle, density and speed of sound at saturation, while d has been deter- 
mined by demanding that the line tension of the semi-infinite system equals 
-^DMC- The conclusions that we draw from the subsequent analysis are: 
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• The energy per particle of two-dimensional slabs approaches that of 
the bulk system slowly, in a linear way. In this limit the radius of 
the slabs diverges as expected. Other physical parameters approach 
the corresponding bulk values much faster. For instance, the chemical 
potential or the thickness have attained their bulk values for inverse 
coverages 5; ~ 1 A^^ 

• For clusters, the bulk limit is harder to attain, and both the energy 
and the chemical potential approach linearly the corresponding limiting 
values. An analysis of the energy per particle in the form of a mass 
formula gives a positive value for the curvature term £c, as one expects 
from the lack of binding energy of the atoms on the surface. 

• Some features of the density profiles of finite clusters are similar to 
their three-dimensional counterparts. For example, the central density 
presents a 'leptodermous' behavior, approaching the bulk saturation 
density from above. 

• Finally, an analysis of these profiles in terms of hydrodynamical equi- 
librium has shown that near the surface of the puddles the pressure 
presents a dip down to negative values. This has been interpreted as 
the effect of the line tension. Hydrodynamic equilibrium demands de 
presence of a force (due to the He-He interaction) that close to the 
surface is proportional to the distance to it. 

Further studies to be performed on *^He systems are the calculation of the 
natural orbits and the subsequent evaluation of the condensate fraction. 
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Appendix A 

Summing up over fermionic 
Matsubara frequencies 



In this appendix we show how to evaluate a summation over Matsubara 
frequencies. First of all, let us shortly introduce the Matsubara technique for 
finite-temperature Green's functions (see [Mat92j). 

At zero temperature, the solution of the Schrodinger equation 



d 



(A.l) 



can be reexpressed in terms of the (zero temperature) Green's function (or 
propagator) 



^0 



(A.2) 



where T is the time ordering operator, 

T[a{t)b{t')] = a{t)b{t')e{t - t') ± b{t')d{t)Q{t' - t) , 

and |\l/o) is the exact ground state of the sytem. The operators in (|A.2|I are 
defined int the Heisenberg picture as 

d{t) ■= e'^^de-'^^ . 

Analogously, for a system at T > 0, we define a propagator as a thermal 
average of the same operator as before: 



G^{x,t; x',t' 



I T 



^{x,t)^\x\t'^ 



(A.3) 
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where the angular brackets stand for the thermal average with the statistical 
operator — exp[— — /xiV)], 



{O) 



Tr [p^O] 



Z = Tr[pT . 



However, the statistical operator satisfies Bloch's equation, namely 

-| = (^ 

which is formaly identical to Schrodinger's equation with the equivalences 

it/h^(3 H^H-pN. 

This suggests to define the finite-temperature Green's function just as the 
zero-temperature one, but performing these replacements everywhere, that 
is, 

(A.4) 



G{xt,x't') := - {T \^{xt)^I}\x't') 

with T — it the so-called "imaginary time". It can be shown that this propa- 
gator (and not the defined above) can be evaluated in a series expansion 
as it can be done with the zero-temperature Green's function. 

We turn now to show three general properties of this G that will be of 
interest. Dropping for simplicity the x-dependence of the propagator, it is 
easy to show that G(r, r') = G{t' — r): 



G'(r,r') = --Tr 



-/SKKt 



-Kr/h Kr 



m^^Ki,r'-r)/n^'^^-Ki,r'-r)/n^ = (7(0, t' - t) , 



where we introduced the grand-canonical Hamiltonian K — H — /iN, and 
used the cyclic property of the trace. 

We can also show that, for fermion, G{—j3h) = — G(0): 



G{-Ph) = --Tr 



-UK 



±G(0) , 



■Ixr 
Z 



where we used again the cyclic property of the trace. 



157 



With these two properties, it is easy to show that G{t) is periodic with 
period ph, and it is useful to introduce its components in Fourier space, 

G{T) = j^J2'"^"^^^''r.). (A.5) 

n 

Then, we have 

n n 

From the condition G{—[3h) = — G(0) we therefore get 

g.;^„/3 ^ _^ ^ ^ + 1)-^ , (A.6) 

which are called fermionic Matsubara frequencies. 

Let us return to the evaluation of a summation over such Matsubara 
frequencies a;„ as the one we encountered in Sect. I1.2T31 evaluation of 

S = y( -T-^-TT - -T-^—fA ■ (A.7) 



ihuJn — Ek ihwn + E, 

The way to perform this kind of sumations is to convert them by means 
of Cauchy's theorem on integrals on the complex plane |Arf95| : 

±,£dzf{z) = ^Res{f{zn)} , (A.8) 

where f{z) is an analytic function inside the region enclosed (in the mathe- 
matically positive sense) by the contour C (see FigETJ except on the points 
{zn}, that are the poles of /; Res{/(z„)} is the residue of / at z„. 

Looking at Eq. (jA.Tjl . we define a function / on the complex plane by 



so that now 



n 

It turns out that tanh(2;) has poles precisely at the points z = Zn, and 
these poles have unit residues. Therefore, f{z) tanh(z) will have poles at the 
desired positions and with the desired residues, and can write 

S = ^^(fdzfiz)tanhiz), (A.ll) 
2 Zm Jc 
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where the contour C encloses the imaginary axis, as shown by the dashed 
line in Fig. lA.lL but not the poles of / on the real axis. 

If we deform C into C + T (see again Fig. lA.Hl . and realize that the 
integral over F will not contribute at |z| ^ cxd because in that limit f{z) oc 



-2 



we get 

S = — — : (f dzf{z) tanh(z) = — — : (b dzf{z)tanh(z) 



2 27TZ J (ji 2 27T2 J (ji 

2 2m ^ 



2m) ^ Res{tanh(z)} = -/?tanh 1^/?^^ , (A.12) 



where we evaluated the integration over C using again Cauchy's theorem but 
now with a ' — ' sign because C is oriented in the mathematically negative 
sense, and we used the fact that /(z)tanh(z) is analytic on the real axis 
except at z = ±PEk/2. Thus, we obtain 

y ( . ^ r. . ^ ^ 1 = -/3tanh ( P^] , (A.13) 

which is the result used in deriving Eqs. ()1.23|1 and ()2.H1 . 



Figure A.l: Complex 2;-plane with the contours C and C' + F. The crosses 
indicate the poles of tanh(2;) while the circles are located at ±Ek. 




Appendix B 



Deduction of the equations of 
motion of spin-/ condensates 



B.l Introduction 

Consider the interaction between two spin-/ atoms: 

2/ 

V{ri - r2) = 5{ri - J] QePf (B.l) 

F=Q 

(B.2) 

where a^r is the s-wave scattering length in the channel of total spin F, and 
M is the mass of the atoms. For bosons, only even F's are allowed; for 
fermions. only odd F's. From now on we shall consider only bosons. For 
example, for / = 1, 2, we have 

V^g2P2 + goPo = 

V = ^4^4 + 52^2 + 50^0 (/ = 2) 

On the other hand, we also have: 

2/ 

/i-/2-5^A^Pf (B.3) 

To show this, consider the addition of two spins / to a total spin F: 

ifl + f2f = fl + fl + 2/i • /2 = 2/(/ + 1) + 2/i • f2 
{fi + f2f = F{F+l) 
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Therefore, 

F(F + l)-2/(/ + l) 

Jl • J2 — ^ ■ 

To take into account all possible channels {F — 0, ...,2/) we must sum over 
F, projecting each contribution with the corresponding projector Pp, that is 

F=0 

For instance, 

/l-/2 = -2Po + P2 

/i • /2 = -6Po - 3P2 + 4P4 (/ = 2) 

Finally, the normalization condition on the projectors provides the iden- 
tity 

2/ 

J2Pf = 1. (B.6) 
From all these equations we can obtain: 



Spin / = 1: 



Therefore, 



V = ^2^2 + ^oPo 

/l-/2 = -2Po + P2 
Po + P2 = 1 



V/ = ^ \ ^ fl -12 = 00 + C2/1 ■ /2 . 



Spin / = 2: 



^ = 94P4: + 92P2 + QqPo 
/i • /2 = -6P0 - 3P2 + 4P4 
Po + P2 + P4 = 1 



Therefore, 



^ = = \ = PoH = — Jrj2 = Co+CiPo+C2/r/2. 
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Thus, we can write the Hamiltonian in the form 

I . . . . f I . . 

a'=-f a',b,b'=-f 

+ XI ^Ci'0+V'^(Po)aa',66'^&'^& [this term only for / = 2] 



a' ,6,6'=-/ 

- E P{^)i^aiFz)aa4a' + Yl #a (^f )aa'4' , \ (B.7) 
a'=-f a'=-f ) 

where ipa is the anihilation operator for an atom in state |l,a), Notice also 
that we have taken into account the effect of a magnetic field in the lowest 
two orders by 



p{z) — g^BB{z) linear Zeeman effect for B — Bz, 

o2 2 

q — — ; — —Bq quadratic Zeeman effect . 



Here g is the gyromagnetic factor for the atom, and huhfs is the hyperfine 
splitting characteristic of the manifold in which the atoms are. 

B.2 Mean- field dynamical equations for a spin- 
1 condensate 

We make the usual mean-field approximation for the field operators, 

^a(r) = (^a(r)) , (B.8) 

The equation of motion for the field il^m can be derived from the Hamiltonia 
by functional differentiation: 

5{H-iiN) 

in- — 



dt 



2M 



V + [Kxt - n 



+ Co 

-p{z){F^^)^ + q{F^ij)^ . (B.9) 
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In the mean field approximation, one can rewrite the several terms in this 



way: 



'^nair) = n{r) , 

a a 

which is the total density of particles at point r. 



ab ab i=x,y,z 

i 

[i = z] {ni - n_i} 

where we have used the matrix expression of the spin operators, i. e. 
the Pauli matrices for spin f — 1, namely 



1 

7^ 



iV2 



1 

1 
1 



-1 



1 

0-1 






1 







1 



-1 



(B.lOa) 
(B.lOb) 
(B.lOc) 



(/) ■ {f'4^)m' above, we use the Pauli matrices to arrive at 



ab 



\m 



111 = 1] + (ni + no - n-i)ipi 

m = 0] (ni + n^i)^o + '2,^i'ipl%l)_i 
-1] + (n_i + no - ni)^_i 
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Collecting the partial results, we get 



ih 



dipi 
dt 



h^ 



2M 



V + Kxt + con(r) 



+ C2 [(ni + no - n_i) ipi + V'o^* i] 
- pVi + #1 



ih 



dt 



h^ 



2M 



^0 



(B.lla) 



(B.llb) 



di/j-i 
dt 



h^ 



2M 



+ Fext + con(r) 



(B.llc) 



We have separated in different lines the contributions: (1) equal for all com- 
ponents, (2) proportional to C2 and related to spin-exchange collisions, and 
(3) terms due to the presence of a magnetic field. Note that, except for the 
linear Zeeman term (oc p oc 5), the equations for m = 1 and m = —1 are 
invariant under the relabelling 1^—1, indicating our freedom to choose 
the orientation of the quantization axis. More generally, the equations are 
invariant under the time-reversal operator, that transforms 1 — 1 and 
B -B. 

These dynamical equations can be written in the shorter form (I5.7jl by 
introducing effective potentials that, in the case of a non- vanishing magnetic 
field, read 

Vl^ = Vext + Con + C2(±ni + no T ^-1) TP + q 



Vr 



eflf 



Kxt + Con + C2(ni + n_i) . 



It is also possible to derive a set of coupled continuity equations from 
((B.lljl . To this end, we define the number density of the m component as 
usual by 

nmir,t) = \i)„Xr,t)\'^ . 
Applying the operator ihd/dt to uq, we get 



ih 



dnp 
dt 



lh^:^%PQ 



dt 
h^ 



2M 
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Now, defining the quantum current of the m component by 
we can rewrite this result as 

which has the form of a continuity equation with a source term given by 
-{2c2/h) ImfToV'o] e ^, with Tq := 2ipiiplip^x. For the m ^ ±1 components, 
a similar treatment yields 

9n±i , ^ . 2c2 r , 

+ V • j±i = — — ImfTiiV'ii] , 

7^1 ^o^^i • 



B.3 Mean-field dynamical equations for a spin- 
2 condensate 

The starting point is the grand-canonical Hamiltonian that, for / = 2 reads: 

2 ^ 2 ^ 

+ XI ;^CoV'a ^^^a'V'a + XI :;C2lpa'4'a' fab ' fa'b'lpb'i'b 



a'=-2 a',b,b'=-2 

a',6,6'=-2 

2 2 

o'=-2 o'=-2 J 

X (aa'|i|&&') (B.13) 

a,a',b,b' 



Again, we derive the equations of motion of the field operators by a 
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functional derivative of the Hamiltonian (|B.12|I : 



ih 



dt 



2M 



V + [Kxt - /i]^ n 



abc 

-p{z) + . 



(B.14) 



As we are working with / = 2, all spin operators can be represented by the 
following 5x5 Pauli matrices: 



/O 

1 






1 










/ 





1 









V 

/ 2 

10 





\ 










1 






1 

0/ 












1 





\ 










1 
-1 0/ 



-2 / 



(B.15a) 



(B.15b) 



(B.15c) 



With them, we can calculate most of the terms in the equations of motion. 
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For example, performing again a mean-field approximation we get 

/ 1 \ 



1 J| 



Jl 



I 



I 1 



{^2 ri ro r-1 r-2) 





V 10/ 



/ ^2 \ 

^1 
^0 

^^-2 J 



^IpO + -0-2 



Analogously, we get 



{Fy) = -iV'2V'i + ^V't ^V'2 - Y 2 V^o 

+ (V'l - V'-i) + #-1 1 Y - ^-2 1 + t^u^-i , (B.17) 



= 21^21' +iV'ir-iV'-ir-2iv-2i% 

Therefore, the term proportional to C2 becomes 



v 



(^-) (^1^0+ ^-2) ^l{Fy) 



- 



(B.18) 



[m = 2] \ 
[m = 1] 

[m = 0] 

[m = — 1] 
[m = -2] y 
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Now we must calculate the Po-terms. Pq projects a two-particle state 
onto its P = component, Pq = |00)(00|. Therefore, I shall use some Clebsh- 
Gordan coefficients to go from the uncoupled basis of states \fm1fm2) to the 
coupled basis |PM): 

\bc) = ^ \FM){FM\fb, fc) ^= J2 \FM){FM\2b, 2c) 

FM FM 

+1 +2 

= |00)(00|26,2c) + |lM)(lM|26,2c) + ^ |2M)(2M|26, 2c) 

A/=-l M=-2 

Applying the definition of Pq, we have 

Po|6c) = |00)(00|26,2c) = |00)(00|6c) 
iPo)ma,bc ■■= {ma\Po\bc) = (2m,2a|00)(00|26,2c) 

Note that {Po)ma,bc = {Po)am,bc- Thus, it is possible the merge the two terms 
coming from the functional derivative into one single term, that cancels the 
2 factor in front of Ci, as we did in Eq. (fRTl . 

The Clebsch-Gordan coefficients for two spin-2 particles that give a total 
P = are: 

, V I ^ mi = ±2, m2 = =f2 , mi = m2 = 
2mi,2m2 00 = <^ , / ^/ 

1 mi = ±l,m2 = =Fl 

Therefore, 

r I (m,a)&(6,c) G {(±2,^2), (0,0)} 
, . ^1 or (m,a)&(6,c) e {(±1,^1)} 
1^0 Wc S _i (^,a)G{(±2,T2),(0,0)}&(6,c)e{(±l,Tl)} 

[ or (m,a)G{(±l,Tl)}&(&,c)G{(±2,T2),(0,0)} 

For example, for m = 2 the Po-terms are: 

E (^o)2a,6c Mc = 5^ (22, 2a|00) (00|26, 2c)^,V'c 

abc 

E7/.:(22,2a|00)5^(00|26,2c)V'6V'c 

a be 

-^C2E<00|2^'2c)^fe^c 

be 



abe abc 
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Note the common factor that will be the same for all m: 

E' = ^(00|26, 2c)V'6Vc = ^ [2^2^-2 + ipl - 2V'iV'-i] ■ 

be 

The other term can be expressed as: 

a * 

Therefore, if we define 

E = i=E' = ^ [2^2^_2 + il^l - 2^i^_i] , 

we can write the equations of motion for the different spinor components 
thus: 



9i 



(B.19a) 



at 



+ C2n 



+ cinV'* lE + (g - 



(B.19b) 



tn^^ = TLoipo 



dt 



+ C2n 



-{F,} (^1 + + t\l-{Fy) (^1 - 



(B.19c) 



+ C2n 



+ cinV'tE + (g + 



(B.lOd) 



dt 



+ C2n [((F,) + ^(F,)) V-i - 2(F,)z/;_2] 
+ Cin?/;;E + 2(2g + p)^/;_2 , 



(B.19e) 
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where Ho = + Vext ~ A'- + con(r) contains the factors common to all 

equations. Explicit expressions for (Fj) (i = x,y,z) and E have been given 
above. 



Forsi altro canterd con rniglior plectio 



FINIS 
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